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Abstract 


In  this  research  project,  we  analyze  the  mathematical  structure  and  numerical  algorithms 
associated  with  Toeplitz  matrices.  Toeplitz  matrices  arise  in  a  number  of  problems  in 
engineering  and  applied  mathematics.  In  many  such  problems,  the  task  is  to  solve  for  certain 
parameters  of  interest  (  such  as  predictor  polynomial,  reflection  coeffients,  and  solution  to  a 
linear  system  )  in  a  computationally  efficient  manner.  Also,  the  numerical  stability  aspects 
of  the  various  algorithms  must  be  examined  from  the  standpoint  of  implementation  using 
finite  precision  arithmetic.  We  have  derived  fast  (  order-recursive  )  and  superfast  (  fast 
Fourier  transform  based  )  algorithms  for  solving  a  Toeplitz  linear  system.  The  algorithms 
reported  here  ate  some  of  the  most  computationally  efficient  algorithms.  Also,  the  numerical 
stability  of  the  split  Levinson  algorithm  is  examined  and  it  is  established  that  it  is  weakly 
stable.  Furthermore  the  various  classical  and  split  Levinson  algorithms  are  studied  for  the 
effects  of  finite  precision  arithmetic.  An  interesting  relationship  between  Levinson  algorithm 

and  stability  tests  for  discrete  systems  is  exploited  to  derive  a  new  computationally  efficient 

/. 


algorithm  for  testing  the  wide  sense  stability  of  discrete  time  systems. 


Chapter  1 


INTRODUCTION 


Toeplitz  matrices  occur  in  a  large  number  of  problems  in  engineering  and  applied  mathe¬ 
matics.  In  problems,  we  have  studied  in  pattern  recognition  and  prediction  and  filtering 
theory,  they  arise  as  the  autocovariance  matrix  of  a  wide  sense  stationary  stochastic  process. 
Some  of  the  important  computational  problems  associated  with  a  Toeplitz  matrix  that  we 
have  studied  are  (  i  )  computation  of  the  associated  predictor  polynomial  and  reflection 
coefficients,  and  (  ii  )  computation  of  the  solution  to  a  Toeplitz  linear  system.  Emphasis 
is  placed  upon  exploiting  the  rich  mathematical  structure  of  the  Toeplitz  matrix  to  derive 
algorithms  having  as  low  arithmetic  complexity  as  possible. 

Another  important  aspect  of  numerical  algorithms  is  their  stability  when  implemented 
using  finite  precision  (  either  fixed  point  or  floating  point  )  arithmetic.  Such  an  analysis  is 
critical  when  the  data  is  perturbed  by  errors  and  the  underlying  problem  is  ill-conditioned. 
The  numerical  stability  of  an  algorithm  determines  its  applicability  for  solving  practical 
problems  of  interest. 
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In  this  research  project,  we  have  focused  on  computational  aspects  of  linear  algebra 
problems  associated  with  Toeplitz  matrices.  The  matrix  structures  considered  are  one  di- 
mensional  as  well  as  two  dimentional.  We  have  derived  fast  (  order-recursive  )  as  well  as 
superfast  (  fast  Fourier  transform  based  )  algorithms  for  solving  a  Toeplitz  linear  system. 
The  algorithms  reported  here  are  some  of  the  most  computationally  efficient  algorithms  for 
solving  Toeplitz  systems.  Also,  the  numerical  stability  of  the  split  Levinson  algorithm  is 
examined  and  it  is  established  that  it  is  weakly  stable.  Furthermore,  the  various  classical 
and  split  algorithms  are  studied  for  the  effects  of  finite  precision  arithmetic.  An  interest¬ 
ing  relatoinship  between  the  Levinson  algorithm  and  stability  tests  for  discrete  systems  is 
exploited  to  derive  a  new  computationally  efficient  algorithm  for  testing  the  wiae  sense 
stability  of  discrete  time  systems. 

The  main  contributions  of  this  research  project  are  as  follows: 

1.  A  new  algorithm  for  testing  the  wide  sense  stability  of  discrete  time  systems  is  derived. 
A  system  is  term  stable  in  the  wide  sense  if  none  of  its  poles  are  outside  the  unit  circle. 
The  new  algorithm  requires  0.25n2  multiplications  (  MULT  )  and  0.5n2  additions  ( 
ADD  ).  This  represents  a  reduction  by  50%  in  the  multiplicative  complexity  over 
previous  algorithms. 

2.  The  numerical  stability  of  the  split  Levinson  algorithm  is  studied  and  various  bounds 
on  the  norm  of  the  residual  vector  are  derived.  It  leads  to  an  important  result  -  the 
split  Levinson  algorithm  is  weakly  stable. 

3.  The  block  Toeplitz  matrices  are  studied  and  new  fast  and  superfast  computationally 
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efficient  algorithms  are  derived  for  solving  a  block  Toeplitz  linear  system.  Additional 
computational  savings  are  realized  when  the  block  Toeplitz  matrix  has  additional 
structure. 

4.  The  effects  of  finite  precision  arithmetic  on  the  computation  of  reflection  coefficients 
using  the  classical  as  well  as  the  new  split  algorithms,  are  studied.  Various  bounds 
are  derived  which  demonstrate  the  propagation  of  errors  from  one  stage  to  the  next 
in  the  order  recursive  algorithm. 

5.  A  new  form  of  three-term  recurrence  relation  is  derived  and  computationally  efficient 
alternatives  to  the  Levinson,  Schur,  lattice  and  normalized  lattice  algorithms  are 
obtained.  We  also  describe  a  new  lattice  realization  of  a  digital  filter  and  a  new 
computationally  efficient  algorithm  for  solving  the  Toeplitz  linear  system. 

This  report  is  organized  as  follows.  In  Chapter  2,  the  Levinson  algorithm  and  its  salient 
feature  are  described.  Also,  based  on  the  Levinson  algorithm,  new  tests  are  derived  for 
testing  the  strict  sense  and  wide  sense  stability  of  a  discrete  time  system.  In  Chapter  3,  the 
numerical  stability  of  the  split  Levinson  algorithm  is  studied.  Chapter  4  deals  with  fast  ( 
order-recursive  )  and  superfast  (  fast  Fourier  transform  based  )  algorithms  for  solving  block 
Toeplitz  linear  systems.  Chapter  5  deals  with  the  finite  precision  analysis  of  classical  as  well 
as  split  algorithms  for  computing  the  reflection  coefficients.  In  Chapter  6,  new  algorithms 
are  derived  for  computing  the  reflection  coeffients  and  subsequently  solving  the  Toeplitz 
linear  system.  Chapter  7  contains  some  conclusions  for  this  research  work. 
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Chapter  2 


On  Fast  Algorithm  For  Testing  The 
Stability  Of  Discrete  Time  Systems 


Stability  is  a  basic  consideration  in  the  design  of  any  system.  Usually  ,  bounded  output  is 
required  for  bounded  input.  The  most  popular  technique  to  test  stability  of  a  given  system  is 
to  locate  the  root  distribution  of  the  characteristic  function  of  the  system.  Several  methods 
for  the  popuse  such  as  Routh  and  Hurwitz  stability  criteria  have  been  proposed.  These 
methods,  however,  are  complicated  in  computation  because  of  their  generality. 

In  this  chapter,  a  fast  and  efficient  algorithm  for  the  testing  stability  will  be  studied.  The 
main  point  is  that  instead  of  computing  the  roots,  a  set  of  parameter  called  the  reflection 
cofficients  are  computed.  A  close  relation  between  the  reflection  cofficients  and  roots  of  a 
given  characteristic  polynomial  is  that  the  roots  are  located  inside  the  unit  circle  if  and 
only  if  the  magnitude  of  each  of  the  reflection  cofficients  is  less  than  one.  The  methods  axe 
based  on  the  Levinson  algorithm  which  is  widely  used  in  digital  signal  processing. 
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2.1  The  Levinson  Algorithm 


The  minimum  least  square  method  is  used  in  the  design  of  discrete  systems  in  applications 
such  as  power  spectrum  estimation, and  speech  prediction.  The  linear  systems  obtained 
often  have  a  special  structure.  This  special  system  is  called  the  Toeplitz  system.  There  are 
efficient  algorithms  for  solving  the  Toeplitz  linear  system  because  of  its  unique  mathemat¬ 
ical  properties.  Levinson  algorithm,  an  efficient  and  fast  method  for  solving  the  Toeplitz 
systems,  was  introduced  by  N.  Levinson  [1].  Toeplitz  systems  occur  widely  in  modern  signal 
processing. 

Suppose  that  a  stationary,  finite  energy  signal  sequence  u(l)  is  given.  If  we  want  to 
predict  the  u(l)  by  its  previous  n  data,  we  form 

«(0  =  £-«••«(*- 0  (2-i) 

«=i 

where  a{,  i  —  1,2,  •  •  •  ,n  are  called  the  prediction  coefficients  and  are  determined  by  mini¬ 
mizing  the  square  of  prediction  error 


<*„  =  (u(f)  -f  ^a,u(/  -  i ) 

\  i=i  / 


(2.2) 


which  results  in  the  following  linear  system  of  equations, 


Tn_i 


/  \ 

(  \ 

Oi 

1 

II 

Cl 

< a" ) 

^  ) 

(2.3) 
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where 


el  =  (1  Cl  c2  •  •  •  Cn)  (2.6) 

and 

sL  =  (lai  •••  an)  (2.7) 

then,  the  minimum  error  an  is  given  by 

an  =  1  +  aci!  +  c2a2  + - \-cnan=£an  (2.8) 

Combining  (2.3)  and  (2.8),  the  normal  equation  can  be  written  as 

Tnan  =  (an0  ...  0)‘  (2.9) 

For  a  given  normal  equation,  one  must  find  the  inverse  of  the  Toeplitz  matrix  to  solve  for 
the  unknown  vector  on.  The  classical  algorithms  for  solving  the  linear  system  axe  Gaussian 
elemination  and  Cholesky  decomposition.  Both  of  them,  however,  are  designed  for  an 
arbitrary  matrix.  The  special  structure  of  Toeplitz  matrix  is  not  exploited  by  these  general 
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methods,  that  is,  even  though  these  classical  methods  can  be  used  to  solve  the  Toeplitz 
system  in  (2.9),  they  are  not  efficient.  Generally  speaking,  the  operation  time  required 
by  Gaussian  elemination  and  Cholesky  decomposition  is  proportional  to  n3.  The  efficient 
algorithm  which  exploits  the  structure  of  the  Toeplitz  matrix  should  have  smaller  order. 

One  of  the  efficient  algorithm  for  solving  Toeplitz  problem  is  Levinson  algorithm.  It  is 
widely  used  for  its  simple  and  efficient  structure.  The  general  Levinson  algorithm  will  be 
studied  in  Chapter  3.  The  details  are  omited  here. 

The  Levinson  algorithm  has  a  form 


PkOCk-l 

J 

II 

(2.10) 

ak,i 

»=0 

=  ak-l,i  +  pkOk-l,k-i 

(2.11) 

Ojfc 

(2.12) 

with  initial  conditions  a^o  =  1,  a k,k  =  Pk  and  ao  =  1- 

In  this  report,  the  polynomial  representation  of  a  sequense  is  also  used.  For  a  given 
vector  a*  =  (a*i0  •  •  •  a*,*)*,  the  polynomial  form  is  defined  as 

k 

a*(*)  =  2  (2.13) 

i=0 

The  Levinson  algorithm  in  polynomial  form  is 

ak{z)  =  ak-i(z)  +  pkz~lak-\{z)  (2.14) 

where  a*_i(z )  denotes  the  reciprocal  polynomial  of  dk-i(z)  and  is  defined  as 

ak(z)  =  z~kak(z~l)  (2.15) 
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After  getting  the  an(z)  by  Levinson  algorithm,  the  inverse  matrix  of  Toeplitz  matrix  Tn 
can  be  obtained  through  the  Gohberg-Semencul  formula 


1  0 


T?  =  a; 


_!  a"'X  1 


0  1  an,i 


0  0  1 


<*n,i  1  0  0 


0  0 


an,n  0 


[  <*n,  1  •  •  •  <*n,n  0  J  [  0  0 

The  general  formula  will  be  derived  in  next  chapter. 

Let  us  define  a  new  matrix  A„  as 


An  = 


1  <*1,1  <*2,2  *  •  *  <*n,n 

0  1  <*2,1  *  ’  -  <*n,n— 1 


0  0  0 


In  the  following  we  show  that 


T-1  =  An 


1  0 


0  c*i 


0  0 


<*n,n— 1 


0  0  <*n,n  ‘  ‘  ‘  <*n,l 

0  0  0  •  •  •  a„,2 


(2.16) 


(2.17) 


(2.18) 
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Proof:  Since  we  have 


TnAn  — 


'  V 
0 
0 


Tn 


/  \ 

01,1 

1 

0 


using  (2.9),  we  have 


v  \° )  V  0  / 


1  0  0 


TnAn  = 


*  at  0 


*  *  OC2 


*  *  * 


1  0  0  •••  0 


a 


"/ 
\  ( 


*  1  0 


*  *  * 


\\ 


‘rijn 


On,n— 1 


On,l 


0  0  0 


)) 


1  0  0  •••  0 


0  10- 


&n  J 


where  the  stars  represent  the  terms  that  we  do  not  care  about.  Multipling  both  side  by 


Anx  we  have 


Tn  =  RnDnA-1 

since  Tn  is  symmetric,  we  have 

r*  =  a-'d^  =  r„ 
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comparing  the  results,  we  see  that  Rn  =  A~ *.  The  proof  is  complete. 

This  result  gets  the  Cholesky  decompositon  of  inverse  Toeplitz  matrix.  The  results 
(2.16)  and  (2.17)  are  different.  The  matrices  in  (2.16)  are  also  Toeplitz.  Based  on  the 
structure,  the  displacement  ranks  of  matrices  which  can  be  represented  by  a  sum  of  two 
product  Toeplitz  matrix  is  studied  by  [2]. 

2.2  The  Classical  Wide  Sense  and  Strict  Sense  Sta¬ 
bility 

The  Levinson  algorithm  introduced  in  last  section  can  be  used  in  application  such  as  linear 
predictive  coding  and  spectral  analysis.  It  is  also  used  in  design  of  ARMA  filter  for  testing 
stability  of  the  filter. 

Stability  is  a  basic  consideration  in  the  design  of  any  system.  For  a  discrete  time  system, 
one  needs  to  check  the  location  of  the  roots  of  the  charecteristic  polynomial  with  repect 
to  the  unit  circle  (  UC  )  for  testing  the  stability  of  the  system  [22,  23].  The  concept  of 
stability  has  been  extended  in  [3].  The  system  is  said  to  be  stable  in  the  strict  sense,  if  all 
the  roots  are  inside  the  UC.  Strict  sense  stability  (SSS)  is  same  as  the  classical  definition 
of  stability  and  can  be  found  in  any  standard  text  book  on  signal  processing.  The  readers 
are  referred  to  [22,  23].  The  system  is  said  to  be  stable  in  the  wide  sense,  if  no  roots  are 
outside  the  UC  (  all  the  roots  lie  inside  or  on  the  UC  ).  The  wide  sense  stability  (  WSS  ) 
is  a  new  idea.  It  was  proposed  after  a  study  of  the  spectral  analysis  of  A  TIMA  filters  and 
adaptive  filters  [4,  24].  In  the  following,  we  will  discuss  <  i  algorithm  for  testing  stability.  It 
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is  based  on  the  Levinson  algorithm.  After  that  we  will  present  a  new  algorithm  for  testing 
the  WSS. 


Let  the  characteristic  polynomial  of  a  given  filter  be 

n 

an(z)  = 

i=0 


(2.19) 


it  is  assumed,  without  any  loss  of  generality,  that  the  first  coefficient  of  an(z)  is  unity. 
From  the  previous  chapter,  the  Levinson  algorithm  can  be  written  as 


ak(z)  =  ak-i(z)  +  pkz  lak-i(z) 


pk  =  ak,k 


(2.20) 

(2.21) 


The  reciprocal  polynomial  of  ak(z)  is 


ak(z)  =  2-1a*-i(z)  +  pkak-i(z ) 


(2.22) 


Solving  (2.20)  and  (2.22)  for  ak-i(z)  leads  to 

.  /  *k(z)  -  Pk&k(z)  ,nnn N 

ak-i[z)  = - - - -  (2.23) 

Equations  (2.21)  and  (2.23)  constitute  the  inverse  Levinson  algorithm.  For  a  given  an(z) 
,  the  reflection  coefficients  pkl  k  =  n,  •  •  • ,  1,  can  be  computed  using  (2.21)  and  (2.23).  It 
has  been  proved  that  the  roots  of  an(z )  are  located  inside  the  UC,  if  and  only  if  (iff)  the 
reflection  coefficients  satisfy 


|pfe|  <  1,  fc  =  l,2,--*,n  (2.24) 

So  for  testing  the  strict  sense  stability  of  a  system,  we  first  compute  the  reflection  coefficients 
using  the  inverse  Levinson  algorithm.  If  the  system  is  stable,  the  condition  (2.24)  should 
be  satisfied.  This  procedure  is  also  called  the  stepdown  procedure  [25]. 
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It  is  clear  from  (2.23)  that  if  \pk\  =  1,  or  if  some  roots  are  located  on  the  unit  circle,  we 
can  not  use  the  classical  stepdown  procedure  to  compute  afc_i(z).  For  the  SSS,  we  claim 
that  the  system  is  unstable  in  this  situation. 

It  is  established  in  [3]  that  if  pk  —  1,  then  ak(z)  must  satisfy 


«*(*)  =  ak(z) 


(2.25) 


and  if  pk  =  —  1,  then  ak(z)  must  satisfy 

ak{z)  =  —ak(z) 


(2.26) 


for  wide  sense  stability.  The  polynomials  which  satisfy  (2.25)  axe  called  symmetric.  The 
polynomials  which  satisfy  (2.26)  are  called  anti-symmetric. 

If  ak(z )  is  symmetric  or  anti-symmetric,  then  it  is  proved  in  [3]  that  afc_!(z)  can  be 
expressed  as 


ak-i(z) 


1  dak(z) 
k  dz~x 


(2.27) 


or  computed  as 

1  *-1 

ak-i(z)  =  1  +  r  ^2(k  -  i)ak<iz~ ‘  (2.28) 

K  t=i 


At  this  point,  we  have  overcome  the  problem  caused  by  =  1  in  classical  step-down 
procedure. 

This  result  leads  to  a  test  for  WSS.  It  can  be  summarized  as:“  the  necessary  and  sufficient 
condition  for  WSS  axe  \pk\  <  1,  k  —  1, 2,  •  •  •  n  (  necessary  condition  for  WSS  )  and  if  =  1, 
then  the  corresponding  polynomial  ak(z)  must  be  either  symmetric  or  anti-symmetric.” 
We  should  note  that  the  symmetric  or  anti-symmetric  polynomials  axe  specified  in  term 
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of  approximate  by  one-half  of  their  coefficients.  In  the  following  section,  we  discuss  the 
tests  for  SSS  and  WSS  based  on  the  recursive  computation  of  symmetric  or  anti-symmetric 
polynomials. 


2.3  New  Step-down  Procedure 

In  the  last  section,  we  discussed  the  classical  Levinson  algorithm  for  testing  the  WSS 
and  SSS.  A  new  fast  algorithm  is  introduced  in  this  section.  This  algorithm  is  called  the 
step-down  precedure.  It  is  based  on  the  recursive  computation  of  either  symmetric  or 
anti-symmetric  polynomials.  We  discuss  the  symmetric  case  first. 

2.3.1  The  Symmetric  Case 

The  new  step-down  procedure  is  based  on  the  recursive  computation  of  the  symmetric 
polynomial  sk(z)  defined  as 

«*(*)  =  a*-i(*)  -(-  z-1afc_i(z)  (2.29) 

It  is  clear  from  (2.29)  that  sk<0  =  3k, k  =  a*-i,o  =  1,  k  —  1,2,  •  •  •  n.  From  (2.23),  we  have 

A  /.a  z~(k~l)[ak{z~l)  -  pkak{z~x)) 

a‘"w  -  - i - 

_  zM*)  -  pkdk(z)] 

1  -Pi 

Dividing  both  sides  by  z  and  adding  the  result  with  (2.23),  we  have 

=  HihliiiM  (2.30) 

1  +  Pk 
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Now  we  derive  the  new  step-down  procedure  [26,  5].  Using  (2.30),  we  can  write 


(1  +  Pk)sk-i(z)  =  ak-i(z)  +  ak-i(z) 

Using  this  equation  with  (2.30),  we  have 

/x  sk(z)  -  ^-1sfc_i(«)(l  +  pk-l) 
ak-i(z)  = - - - 1-z— - 

and 

,  ,  x  Sfc-l(g)(l  +  pk- l)  -  sk(z) 


Substituting  them  into  (2.23)  and  rearranging  the  terms,  we  have 

Sfc-i(z)  =  alXz[{  1  +  z~1)sk(z)  -  sjfc+i(^)] 


(2.31) 


where 


=  1  +  8jt,i  —  Sfc+ l,l 


(2.32) 


and 

Pk-,  =  -1  +  r— ^ —  (2-33) 

In  the  above,  (2.32)  is  obtained  by  noting  that  sk-i,o  =  1.  Given  the  polynomial  an(z), 
the  initial  polynomials  Sfc+i(z)  and  sn(x)  axe  obtained  from  (2.29)  and  (2.30)  by  setting 
k  —  n  +  1  and  k  =  n  respectively,,  where  pn  is  obtained  from  pn  =  an,„. 


2.3.2  The  Anti-symmetric  Case 

In  this  case,  the  new  step-down  procedure  is  based  on  the  recursive  computation  of  the 
anti-symmetric  polynomial  tk(z)  defined  as 

tk(z)  =  o*_i(x)  -  z-'ak-iiz)  (2.34) 
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Once  again,  we  have  tkt0  =  —tk,k  =  1?  &  =  1,2,  •  ■  •  n.  The  derivation  is  very  similar  to  the 
symmetric  case  and  the  various  expressions  can  be  summarized  as  follows 

M»)  =  ’  (2-35) 

tk-i(z)  =  Pklz[(l  +  z~')tk(z)  -  tfc+i(z)]  (2.36) 

where 

Pk  —  1  +  tk<  i  —  tk+ i,i  (2.37) 

and 

/>*-!  =  1  -  rr-  (2-38) 

1  +  pk 

Given  the  polynomial  an(z),  the  initial  polynomials  tn+i(z)  and  tn(z)  are  obtained  from 
(2.34)  and  (2.35)  by  setting  k  =  n  +  1  and  k  =  n  respectively,  where  pn  is  obtained  from 
Pn  ~  an,n-  Since  either  one  of  (2.31)-(2.33)  or  (2.35)-(2.38)  may  be  employed  for  the 
reflection  coefficient  computation,  they  may  also  be  used  as  the  test  for  SSS.  The  reflection 
coefficients  are  computed  using  either  (2.31)-(2.33)  or  (2.35)-(2.38)  and  if  |/?fc|  <  1,  k  = 
1, 2,  •  •  •  n,  then  the  system  is  stable  in  the  strict  sense. 

We  will  give  an  example  to  describe  the  procedure  of  the  algorithms.  We  use  the 
symmetric  case,  the  anti-symmetric  case  is  almost  the  same  as  symmetric  case. 

Example  1.  Let  Zn  denote  the  vector  Zn  =  (1  z~x  •••  z~n)*.  Consider  the  transfer 
function  aA(z)  =  (1  1.6  0.11  —  0.844  —  0.336)Z4.  Here,  n  =  4,  and  p4  =  —0.336.  We  use 
(2.39)  and  (2.40)  to  compute  ss(z)  as,  s$(z)  =  (1  1.264  —  0.734  —  0.734  1.264  1)Z5  and 
s4(z)  =  (1  1.1386  0.3313  1.386  1)Z4.  Employing  the  recurrence  in  (2.31)-(2.33),  we  obtain 

a4  =  0.8746 
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0.3454 


pz  — 

s3(z)  =  (1  2.52  2.52  1  )Z3 
a3  =  2.3814 
p2  =  0.7701 
s2{z)  =  (1  1.9773  1  )Z2 
a2  =  0.4573 
px  =  0.9890 
Sl(z)  =  (11  )Z 

Since,  \pk\  <  1,  k  =  1,2, 3,4,  the  polynomial  is  stable  in  the  strict  sense.  The  roots  of  a4(z) 
are  z\  =  0.2,  z2  =  —0.3,  z3  =  —0.7,  z4  =  —0.8.  They  are  all  inside  the  unit  circle. 

We  observe  that  if  one  of  the  reflection  coefficients,  say  />*  is  such  that  |p*|  =  1,  then 
the  step-down  procedures  can  not  be  used  for  the  reflection  coefficient  computation.  In  the 
next  section,  we  derive  alternative  techniques  to  incorporate  such  cases  referred  to  as  the 
singular  cases. 

2.3.3  The  Singular  Cases 

The  condition  |p*|  <  1  is  only  a  necessary  but  not  a  sufficient  condition  for  WSS.  In 
addition,  the  new  step-down  procedures  as  outlined  in  last  section  can  not  used  if  \pk\  =  1. 
In  this  section,  we  derive  modifications  to  the  new  step-down  procedures  to  incorporate  the 
case  \pk\  —  1.  This  implies  that  either  pk  =  1  or  pk  =  —  1.  In  the  following,  we  discuss  the 
cases  pt,  =  —  1  and  pk  =  1  separately.  We  begin,  our  analysis  with  the  case  pk  =  —1. 
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i).  pk  =  -1 

In  this  case,  for  testing  the  WSS,  we  need  to  check  if  a^(z )  is  anti-symmetric  or  not. 
The  following  lemmas  play  a  key  role  for  such  a  check. 


Lemma  1 


Whenever  pk  =  —  1  and  |/?jt+1|  ^  1,  the  polynomial  ak{z)  is  anti-symmetric  iff 


(1  +  2  1)sfc+i(2)  -  Sk+ffz)  =  0 


(2.39) 


Proof:  Iterating  (2.31)  one  step  forward,  we  have 

Sfc(z)  =  a^z^l  +  2_1)sfc+1(z)  -  Sfc+2(z)] 

subsituting  Sfc(z),  in  terms  of  ak(z)  and  ak(z),  using  (2.28),  we  obtain 

ak(z)  +  ak(z)  =  Z(  1  -  /»fc+l)-1[(l  +  Z~1)sk+l  ~  5fc+2(2)] 

Using  the  above  expression,  the  statement  of  the  lemma  follows  in  a  straightforward  manner. 
Also  if  ak(z)  =  — dfc(z),  then  using  (2.29),  we  get 

(2-40) 

Lemma  2  Whenever  pk  =.— 1  and  |^it+i|  ^  1,  the  polynomial  ak{z)  is  anti-symmetric  iff 

(l  +  z-1)<fc(z)-^+1(z)  =  0  (2.41) 


Proof:  Using  (2.34)  and  (2.35),  we  have 


*fc+i(2) 

tk(z) 


ak{z)  ~  z  lak{z) 

ak(z)  -  ak(z) 

1  -  Pk 
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Solving  for  ak(z)  and  ak(z),  we  have 


ak(z)  =  (1-2  *)  l[tk+1(z)  -  ztk(z)(l  -  pk)\ 
ak(z)  =  (1  -  2"1)‘1[ifc+1(2)  -  tk(z){l  -  pk)} 

Adding  them  and  setting  pk  =  —1,  we  have 

ak{z)  +  ak(z)  =  -2(1  -  z-1)-1[(l  +  z~x)tk(z)  -  4+i(2)] 

and  the  lemma  follows. 

Also,  if  ak(z)  =  —  ak(z),  then  using  (2.35),  we  obtain 

ak(z)  =  tk(z)  (2.42) 

ii).  pk  =  1 

In  this  case,  for  testing  the  WSS,  we  need  to  check  if  ak(z )  is  symmetric  or  not.  The 
following  lemmas  play  a  key  role  for  such  a  check. 

Lemma  3  Whenever  pk  =  1  and  |/9jt+i|  1,  the  polynomial  ak(z)  is  symmetric  iff, 

(1  +  z~x)sk{z)  -  sk+i(z)  =  0  (2.43) 

If  (2.43)  holds,  then  we  also  have 

ak(z)  =  sk(z)  (2.44) 

Lemma  4  Whenever  pk  =  1  and  |/>*+i|  ^  1,  the  polynomial  ak(z)  is  symmetric  iff, 

(1  +  z_1)*fc+1(2)  -  <*+2(*)  =  0  (2.45) 
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If  (2.45)  holds,  then  we  also  have 


ak{z)  =  {l-  z  ')  1tk+i(z)  (2.46) 

The  proofs  for  Lemma  3  and  Lemma  4  are  similar  to  the  proofs  for  Lemma  2  and  Lemma  1 
respectively.  On  the  basis  of  the  above  analysis  and  the  test  for  WSS  using  dk{z),  the  tests 
for  WSS  that  employ  the  polynomials  5^(2)  or  tk(z )  can  be  described  in  a  straightforward 
manner.  In  the  following,  we  summarize  the  test  for  WSS  using  sk{z).  The  test  for  WSS 
using  tk{z)  may  be  described  in  a  similar  manner. 

Test  for  W55  using  sk(z) 

Step  1.  Use  the  stepdown  procedure  given  in  (2.31)-(2.33)  and  compute -p*,  k  =  n,n— 1,  •  •  • 
If  l/Jjtl  >  1,  then  declare  “  unstable  ”  and  stop. 

If  pk  —  —1,  Go  to  step  2.  If  pk  —  1,  Go  to  step  3. 

Step  2.  Check  if  (2.39)  is  satisfied. 

If  not,  then  declare  “  unstable  ”  and  stop. 

If  yes,  compute  a*^.?)  and  ak-i(z)  using  (2.40)  and  (2.28)  respectively.  Go  to  Step  1. 

Step  3.  Check  if  (2.43)  is  satisfied. 

If  not,  then  declare  “  unstable  ”  and  stop. 

If  yes,  then  obtain  ak(z)  and  a*_i(z)  unsing  (2.44)  and  (2.28)  respectively.  Go  to 
Step  1. 
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2.3.4  Computational  Complexity  Consideraton 

Using  the  symmetry  properties  of  the  Sk(z),  it  may  be  established  that  if  none  of  the 
reflection  coefficients  have  unit  magnitude,  then  the  new  step-down  procedure  requires 
0.25n2  +  n  MULT  and  0.5n2  +  4 n  ADD.  If  one  of  the  reflection  coefficients,  say  pk  =  —1, 
then  it  requires  an  additional  2k  MULT  and  0.5 k  ADD.  Similarly,  additional  2k  MULT  and 
0.5 k  ADD  are  required  if  pk  =  1.  Similarly,  using  the  symmetry  properties  of  the  tfc(z),  it 
may  be  established  that  if  none  of  the  reflection  coefficients  have  unit  magnitude,  then  the 
step-down  procedure  requires  0.25n2  +  0.5n  MULT  and  0.5n2  +  3n  ADD.  If  pk  =  —1,  then 
it  requires  an  additional  k  MULT  and  0.5fc  ADD.  If  pk  =  1,  then  it  requires  an  additional 
2k  MULT  and  1.5 A:  ADD.  Thus,  the  procedures  developed  here  require  one-half  the  number 
of  MULT  and  the  same  number  of  ADD  as  the  procedures  described  in  last  section.  The 
computational  savings  in  this  work  stem  basically  from  a  more  efficient  computation  of 
reflection  coefficients.  Also,  the  computational  complexity  of  the  algorithms  presented 
here  is  the  same  as  the  algorithms  presented  in  [27].  However,  the  methodology  and  the 
technique  for  testing  the  stability  in  the  algorithms  presented  here  and  those  in  [27]  are 
entirely  different.  In  the  following,  we  give  severed  examples  to  illustrate  the  algorithm. 

Example  2.  Let  a4(z)  —  (1  0.4  0.48  0.68  —  0.4)Z4.  Here,  n  =  4,  and  p4  =  —0.4.  We  use 
(2.29)  and  (2.30)  to  compute  s$(z)  and  $4(2)  as 

3S(z)  =  (10  1.16  1.16  0  1  )Z5 

s4(z)  =  (1  1.8  1.6  1.8  1)Z4 
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Employing  the  recurrence  in  (2.31)— (2.33),  we  obtain 


£*4  =  2.8 
p3  =  1 

s3(z)  =  (1  0.8  0.8  1  )Z3 


Since  p3  =  1,  we  check  that  Lemma  3  is  satisfied.  Thus  a3(z )  =  s3(z),  and  using  (2.28),  we 
obtain 


a2{z) 


P2 

s2(z) 


(1 


i.e  os  2 

3  3 


0.8 

3 


a  j§ 


2 


Using  (2.31)-(2.33)  once  again,  we  obtain 


a2 


Pi 


19.8 

19 

8_ 
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si(z)  =  (1  1  )Z 


Since  \pk\  <  1,  k  —  1,2, 3,4  and  Lemma  3  is  satisfied  for  pz  =  1,  04(2)  is  stable  in  wide 
sense.  The  roots  of  a4(z)  are  zx  =  0.4, z2  =  -1,  *3,  z4  =  -0.1  ±0.99498;  and  the  magnitudes 
of  z2,  z3  and  z4  are  one. 

Example  3.  Let 

a4{z)  =  (1  0.5  -  1.04  -  0.76  0.3)Z4 
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Here,  n  =  4  and  p4  =  0.3.  We  use  (2.29)  and  (2.30)  to  compute  s5(z)  and  s4{z)  as 


ss(z)  =  (1  0.8  -  1.8  -  1.8  0.8  1)ZS 
aA(z)  =  (1  -0.2  -  1.6  -0.2  1)Z4 


In  this  case,  ot4  =  0  implying  that  p3  =  —1.  We  check  that  Lemma  1  is  satisfied.  Obtain 
a3(z)  using  (2.40)  as 

.  M*)  =  =  C1  °-8  -  °-8  ~  W3 

1  —  z 

and  obtain  a3{z)  using  (2.34)  as 


Finally, 


a2(z)  =  (1  ^ 


4.4 


0*3  = 

3 

0.8 

P2  = 

3 

s2(z)  = 

0  £  U* 

72.2 

0-2  = 

33 

8 

Pi  = 

11 

Si(z)  = 

(1  1  )Z 

Once  again,  the  conditions  for  WSS  are  satisfied.  The  roots  of  a4(z)  axe  z\  —  0.3,  z3  =  1, 
z3,z4  =  0.9  ±  0.43589j.  Note  that  j^l  =  |*3|  =  \z<\  =  1. 

Example  4.  Let 

a5(z)  =  (1  1.3  -  2.6  -  1.9  1.4  0.8)Z5 
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Here,  n  —  5,  and  p$  =  0.8.  We  use  (2.29)  and  (2.30)  to  compute  s6(.z)  and  s5(z)  as, 


s6(z)  =  (12.1  -1.2  -3.8  -1.2  2.11  )Z6 
s5(z )  =  (11.5  -2.5  —  2.5  1.5  1)Z5 

Employing  the  iecurrence  in  (2.31)-(2.33),  we  obtain 

£*5  =  0.4 

p4  =  1 

s4(z)  =  (1  0.5  -3  0.5  1)Z4 

Since  p4  =  1,  we  check  that  Lemma  3  is  satisfied.  Thus,  a4(z )  =  s4(.z)  and  using  (2.34),  we 
obtain 

a3(z)  =  (1  ^  -  1.5  i )Z 3 
1 

*  =  8 

s3(*)  =  (1-1-1  1  )Z3 

Using  (2.31)-(2.33),  once  again,  we  obtain 


<*3 

P2 

s2(z) 


-0.5 

11 

7 

(1  —  2  l)Z2 


Now,  a2  =  0  imply  p\  =  —1.  We  check  that  Lemma  1  is  satisfied,  and  obtain  a\(z )  as 

«.(*)  =  =  (i  - 1 )z 
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Finally,  si(z)  =  (1  1  )Z. 

Since,  \p2\  =  y  >  1,  the  conditions  for  WSS  are  not  satisfied  and  as(z )  is  unstable.  The 
roots  of  a$(z)  are  z\  —  1,  z2  =  1,  23  —  0.5,  z\  —  —2,  z$  —  0.8. 

2,4  Further  analysis  of  the  sigular  cases 

Here,  we  analyze  the  singular  cases  in  terms  of  roots  of  the  polynomial  an(z )  that  lie  on 
unit  circle  or  occur  in  reciprocal  pairs.  Let  u1?  u2,  •••,«!  be  the  roots  on  unit  circle,  that 
is  1^1  =  i,  j  =  1, 2,  •  *  • ,  /  and  n,  r2,  •  •  •,  rm  be  the  roots  occuring  in  reciprocal  pairs.  Also, 
let  ui+m(z)  be  a  polynomial  of  degree  l  +  m  having  m,  u2,  •  •  *,  u/,  and  n,  r2,  •  •  •,  rm  as  its 
roots,  such  that  ul+m<0  =  1.  Clearly,  ul+m(z)  is  a  factor  of  an(z),  and  we  can  write 

an(z)  =  ui+m(z)rn_/_m(z)  (2-47) 

It  is  straightforward  to  verify  that  if  1  is  a  root  of  U[+m(z)  of  even  multiplicity,  then 
is  a  symmetric  polynomial,  and  we  have 

an(z)  =  ui+m(z)rn-i-m(z)  (2.48) 

Similarly,  if  1  is  a  root  of  ui+m(z)  of  odd  multiplicity,  then  u/+m(z)  is  an  anti- symmetric 
polynomial,  and  we  have 

an(z)  =  - ui+m(z)rn-i-m(z )  (2.49) 

In  either  case,  u,+m(z)  is  a  factor  of  on(*)  also,  therby  implying  from  (2.23)  that  u/+m(z) 
is  a  factor  of  an-i(z)  and  so  on.  Therefore,  we  obtain 

ai+m(z )  =  ui+m{z )  (2.50) 


24 


Also,  pi+m  =  1  if  a/+m(z)  is  symmetric  and  pl+m  =  —1  if  aj+m(z)  is  anti-symmetric.  In 
either  case,  =  1-  Consequently,  the  singular  cases  detect  and  check  for  the  presence 

of  roots  on  the  unit  circle  and  roots  occuring  in  reciprocal  pairs. 

Using  (2.47),  (2.48)  and  the  definition  of  sn+i(z)  and  sn(z)  in  (2.29)  and  (2.30),  we  see 
that  tq+m(z)  is  a  factor  of  sn+i(z)  and  sn(z).  Therefore,  (2.31)  implies  that  it  is  also  of 
•Sn-i(^)  and  so  on.  Consequently,  from  (2.29),  we  obtain 

Si+m(z)  =  u/+m(z)  =  ai+m(z)  (2.51) 

if  ai+m(z )  is  symmetric,  or 

s/+m+i(2)  =  (1  -  z-1)tq+m(z) 

=  (l-z->/+m(z)  (2.52) 

if  ai+m(z )  is  anti-symmetric. 

Similar  statements  hold  for  polynomials  f*(z),  k  —  n,n  —  1,  •  •  • ,  /.  It  is  interesting  to 
note  here  that  the  singular  cases  will  occur  as  many  times  as  the  largest  multiplicity  of  a 
root  on  the  unit  circle  or  a  root  occurring  in  reciprocal  pairs  [27].  This  is  due  to  the  fact 
that  a  root  of  a*,  of  multiplicity  gamma  (  7  >  1  )  is  also  a  root  of  of  multiplicity  7  —  1. 

Numerical  properties  ( in  terms  of  stability  and  sensitivity  )  of  the  algorithms  described 
here  as  compared  to  their  classical  forms  axe  an  important  issue  and  needs  further  unves- 
tigation.  However,  such  an  analysis  goes  beyond  the  scope  of  this  work.  Finally,  we  end 
this  section  by  stating  that  computationally  efficient  tests  for  SSS  and  WSS  based  on  the 
stepdown  procedure  may  also  be  described  for  polynomials  with  complex  coefficients.  The 
details  of  such  an  analysis  will  be  presented  elsewhere. 


25 


In  this  chaper,  we  introduced  the  Levinson  algorithm  for  solving  Toeplitz  systems.  We 
also  studied  the  computationally  efficient  procedures  for  testing  the  wide  sense  stability 
and  strict  sense  stability  for  polynomials  with  real  coefficients.  The  procedures  are  based 
on  efficient  procedures  for  computing  the  associated  reflection  coefficients. 


Chapter  3 


Stability  of  Algorithm 


The  algorithm^  for  solving  linear  system  usually  involve  a  lot  of  computations.  The  precision 
associated  with  the  initial  values  of  problem  is  limited,  and  therefore,  that  is,  errors  are 
inevitable.  For  some  problem,  the  algorithms  give  incorrect  results  even  though  the  errors 
in  the  initial  values  are  very  small.  This  problem  has  attracted  a  great  deal  of  attention.  In 
fact,  an  algorithm  is  termed  stable  if  it  gives  correct  results  within  the  required  precision 
range  for  small  errors  in  the  initial  values.  If  an  algorithm  is  stable,  for  small  errors,  one 
can  guarantee  that  the  result  also  has  small  errors. 

In  this  chapter,  we  explore  the  numerical  stability  properties  of  the  recently  reported 
split  Levinson  algorithm  for  the  computing  the  predictor  polynomial  associated  with  a 
positive-definite  real  symmetric  Toeplitz  matrix.  Various  bounds  on  the  residual  vector 
are  derived  for  the  fixed-point  and  floating-point  implementation  of  the  algorithm.  These 
bounds  are  similar  in  form  to  the  bounds  derived  by  Cybenko  for  the  Levinson  algorithm 
and  are  obtained  by  converting  a  three-term  recurrence  for  the  error  vector  to  an  equivalent 
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two-term  recurrence.  The  main  conclusion  of  the  paper  is  that  the  split  Levinson  algorithm 
is  weakly  stable. 

3.1  Stability  of  Levinson  Algorithm 

In  the  chaper  2,  we  have  studied  the  Levinson  algorithm.  It  is  widely  used  in  the  signal 
processing  applications  such  as  speech  processing,  estimation  of  power  spectrum  and  so  on. 
However,  the  stability  properties  of  the  algorithm  are  not  completely  known.  On  one  hand, 
the  theoretical  results  showed  that  the  algorithm  is  sensitive  to  errors  [13]  while  on  the  other 
hand,  the  practical  results  showed  that  the  algorithm  gives  accurate  results  [14].  Recently, 
the  stability  of  Levinson  algorithm  (  it  is  at  least  as  stable  Cholesky  method  )  was  proved. 
This  proof  ended  the  controversy  about  stability  of  Levinson  algorithm.  It  is  shown  that 
the  algorithm  is  ill-conditioned  when  pk,  the  reflection  coefficients,  are  near  the  unit  circle. 
In  this  case,  the  algorithm  is  far  more  sensitive  to  errors.  The  analysis,  strictly  speeking, 
falls  in  the  category  called  the  weak  stability  analysis.  The  concept  of  weak  stability  was 
introduced  in  [7].  According  to  Buch,  the  result  in  [6]  showed  that  Levinson  algorithm  is 
weakly  stable. 

The  proof  of  strictly  stability  of  Levinson  algorithm  is  still  lacking.  Since  errors  of  the 
algorithm  are  intricately  interleaved,  they  axe  very  difficult  to  manipulate.  In  this  section, 
we  will  study  stability  aspects  of  Levinson  algorithm. 

In  Chapter  2,  the  Levinson  algorithm  was  studied.  For  convenience,  we  write  the  algo- 
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rithm  in  the  vector  form  as  follows, 


Pi  = 


g,  = 


OLi  = 


«t-l 

(  fe-i  +  Ai-i  J 

(l  -  Pi)  *  =  1,2. 


',** 


(3.1) 

(3.2) 

(3.3) 


The  initial  condition  is  <*0  =  1.  Here,  pi  is  the  called  reflection  coefficient,  and  the  notation 
represents  the  operation  of  symmetric  reflection  of  a,,  i.e. 

Qu  =  (<*»,«  <*%,«- i  •  •  •  ai) 

The  necessary  and  sufficient  condition  for  Tn  to  be  positive  definite  is  \p,\  <  1.  This  also 
guarantees  that  the  HR  filter,  whose  characteristic  polynomial  is  1  +  £)’=1  g.jz--7,  is  stable. 
The  error  a,-  decreases  monotonically  to  zero  as  i  is  increased. 

It  is  well  known  that  the  condition  of  matrix  A  is  very  important  in  solving  Ax  =  b. 
The  error  of  x  can  be  bounded  as 


\x  —  x 


k(A)\\E\\ 


+ 


*(A)||e|| 


(3.4) 


Hall  "(i-fc(A)|jf|()||Ai|  (i-fc(A)g|)||4|| 

where  x  represents  the  solution  obtained  by  a  computer.  In  this  work,  we  assume  that 
the  matrix  A  is  correct,  i.e.  ||2?||  =  0  and  we  use  the  maximum  column  sum  matrix  norm. 
The  reader  is  referred  to  [15]  for  a  discussion  of  the  norm.  Under  these  assumptions^  3.4) 
becomes 

lla -a II  <  l(a)M  ro  k\ 

Itell  -  (  ’lliil  (  5) 

where  e  is  residual  vector. 
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For  the  Toeplitz  matrix  T„,  Cybenko  has  obtained  the  bound  on  ||Tn  1||  as  [Cybenko, 


1980], 

<36) 

Since  the  |c;|  <  1,  the  norm  of  Tn  satisfies  1  <  j|Tn||  <  n,  and  the  bound  on  condition  of 
Toeplitz  matrix  Tn  is  given  by 


max 


1  _ 1 _ 1  ^  .  ,rp  \  ^  tt  f  ~l~  lQjl 

.«»-i  ’  riy=i  (1  -  Pi)  }  -  n  ”  ”  j5i  1  -  \ai  \ ' 


(3.7) 


This  is  a  very  important  result.  It  reveals  that  Yule- Walker  equation  is  ill-conditioned  for 
|/),  |  near  one  or  for  large  n.  If  this  happens,  one  must  be  very  careful  in  the  computations. 
The  result  is  very  sensitive  to  error.  That  is  why  Box  and  Jenkins  warned  that  the  algorithm 
is  sensitive  to  error  when  n  is  large  [13]. 

Cybenko  has  obtained  bound  on  the  norm  of  residual  for  Levinson  algorithm,  which  are, 


n 

116.11  2  a<t  +  —  +  ">  no  +  M)  +  0(  A2) 


(3.8) 


if  fixed-point  arithmetic  is  used,  and 

116.11  <  A(y  +  lln)  (fid  +  \Pi\)  -  l)  +  0(A2)  (3.9) 

if  floating-point  arithmetic  is  used.  These  results  tell  us  that  Levinson  algorithm  gives 
a  small  residual  if  n  is  not  large  and  A  is  small.  In  other  words,  Levinson  algorithm  is 
weakly  stable  as  pointed  out  by  Bunch  [7].  The  implication  is  that  the  residual  is  small 
even  though  |p,-|  is  near  to  one  if  n  and  A  are  properly  chosen.  As  a  result,  engineers  use 
Levinson  algorithm  in  speech  processing  and  obtain  reasonable  results. 


30 


Markel  and  Gray  have  simulated  the  Levinson  algorithm  in  the  fixed-point  arithmetic 
[14].  Their  results  support  our  argument  here.  Their  results  showed  that  for  minimum 
wordlength,  it  is  necessary  to  use  preemphasis  and  a  sampling  rate  F„  as  low  as  can  be 
accepted.  Preemphsis  can  reduce  |p,|  and  increase  a,-.  The  low  Fs  implies  only  a  few 
sampling  data  points.  They  showed  that  as  Fs  is  increased,  wordlength  must  also  be 
increased.  Recall  that  the  bound  on  quantization  A  is  proportion  to  2~ 0 . 

Another  point  we  note  is  that,  in  practice,  we  require  that  linear  prediction  give  a 
good  result,  that  is,  the  prediction  error  a*  be  small.  However,  smaller  or,-  means  an  ill- 
conditioned  problem,  and  the  computed  result  may  be  more  sensitive  to  quantization  error. 
Consequently,  we  can  not  obtain  an  arbitrarily  small  value  of  a,  and  a  compromise  must 
be  made. 


3.2  Stability  of  Split  Levinson  Algorithm 


Split  Levinson  algorithm(  SLA  )  was  introduced  recently  in  order  to  increase  the  computa¬ 
tion  speed.  It  requires  only  0.5n2  +  0(n)  arithmetic  operations  as  it  manipulates  symmetric 
or  anti-symmetric  vector.  In  the  following,  we  state  some  relevant  results  which  will  be  used 
in  this  report.  The  details  are  given  in  [8].  We  know  that  the  Levinson  algorithm  recursively 
solves  the  system, 


U+l 


=  Ko  0)' 


(3.10) 
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and  the  recurrence  is 


/  ^  .  \ 
fit-  1  +  Pidi- 1 


V 


pi 


) 


Note  that  the  vector  in  this  chapter  is  defined  as 

di  =  {ai,  1  -  ’  ai,i) 

Define  a  new  vector 

dj'-l  =  fii-1  +  fit-1 

After  some  manipulation,  we  get  the  SLA  as, 


-di 


(  \ 

(  \ 

di- 1 

i 

1.  = 

v  1  > 

+ 

l  \ 
1 


Si-2 


\  1  ) 


si  =  2  —  2  di 


= 


/  ,  J  \ 

1  +  S\  —  U2 

1  +  Si  —  di 
i-i 


ti-l 


m=l 


I  "h  ^  Cm^i- l,m  H”  Ct 

di 

Pi  =  1  - 


di 


(3.11) 


(3.12) 


(3.13) 


(3.14) 

(3.15) 

(3.16) 

(3.17) 

(3.18) 


1  +  Pi-1 

with  initial  conditions  t0  —  1,  po  =  0.  Since  |p,|  <  1,  0  <  di  <  4.  Observe  that  s,  is  a 
symmetric  vector.  In  fact,  the  SLA  solves  the  linear  system, 


T’n-lln-l  =  -(cn— 1  +&,-i) 


(3.19) 
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In  this  chapter,  we  study  the  numerical  stability  properties  of  the  SLA.  In  this  section, 
we  state  the  results  and  discuss  the  condition  under  which  the  algorithm  is  well-conditioned. 
We  give  the  various  proofs  in  the  following  sections. 

The  residual  for  the  algorithm  is  defined  as 

Ln  =  Tnan  (3.20) 

The  bound  on  the  residual  for  SLA  is 

IM  <  A  no  +  M  (3.21) 

for  the  fixed-point  arithmetic,  and 

II Ln\\  <  16A  (n3  +  n2  +  n)  JJ  C1  +  bm|)  (3.22) 

m=l 

for  the  floating-point  arithmetic. 

From  these  bounds,  we  see  that,  similar  to  the  Levinson  algorithm,  the  SLA  is  also 
weakly  stable,  if  the  order  n  is  not  large  and  |p,|  is  not  close  to  one.  The  bounds  obtained 
here  are  almost  the  same  as  those  for  Levinson  algorithm  obtained  by  Cybenko.  The  main 
difference  between  SLA  and  Levinson  is  the  power  of  n  in  the  bounds.  Can  we  obtain  a 
bound  which  has  same  power  of  n  as  Levinson  algorithm?  The  answer  is  yes,  but  the  bound 
is  too  loose  and  will  not  be  pursued  here.  We  will  derive  the  bounds  in  (3.21)  and  (3.22) 
in  the  following  sections.  An  important  question  arises  at  this  stage  — ‘  why  the  power  of 
n  is  different  in  the  bounds  for  the  Levinson  algorthm  and  SLA’.  The  reason  is  that  after 
computing  the  symmetric  vector  sn,  we  need  to  go  back  to  prediction  coefficients  an.  This 
,-equires  a  division  by  (1  —  z -1),  if  we  use  polynomial  notation  .  If  we  w;h>i  1  •  >  design  an  AR 
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filter  usi_ig  SLA  (  which  we  can  do  by  using  Levinson  algorithm  ),  the  AR  filter  will  have  a 
terminal  stage  for  division  by  (1  —  z~x).  This  stages  corresponds  to  a  marginally  stable  HR 
filter  of  order  one  having  a  pole  at  z  =  1.  The  divisor  (1  —  z~x)  plays  a  very  important  role 
in  the  SLA.  It  is  the  bridge  linking  the  prediction  coefficients  and  symmetric  vector.  Given 
the  symmetric  vector,  the  only  way  to  get  prediction  coefficients  is  the  bridge.  However, 
since  there  are  errors  in  the  symmetric  vector  due  to  the  finite  wordlength  of  the  computer, 
the  polynomial  formed  by  the  computed  vector  could  not  be  divided  by  1  —  z~x.  The  bridge 
is  broken  in  a  finite  precision  implementaion  of  the  SLA.  Does  it  mean  that  the  SLA  is 
not  stable?  We  point  out  that  for  anti-symmetric  form  of  the  SLA,  the  same  result  can  be 
obtained. 

At  this  point,  we  can  compare  the  Levinson  algorithm  and  SLA.  The  Levinson  algorithm 
can  be  used  to  implement  an  AR  or  ARMA  filter,  while  the  SLA  can  not  be  used.  The 
Levinson  algorithm  gives  directly  results  of  all  the  subsystems  during  the  intermediate 
computations,  the  SLA  does  not.  The  computation  time  of  Levinson  algorithm  is  n2  +  0(n ), 
the  time  of  SLA  is  0.5n2  +  0(n).  However,  from  the  previous  discussion,  we  know  that  the 
problem  is  ill-conditioned  when  n  is  large.  This  means  that  both  the  algorithms  may  be  used 
only  when  n  is  small.  In  such  a  situation,  SLA  may  not  save  computation  time  significantly. 
It  appears  that  Levinson  algorithm  is  safer  than  SLA. 
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3.3  Residual  of  Split  Levinson  Algorithm 


It  is  well  known  that  errors  are  produced  by  the  computer  due  to  finite  wordlengths.  The 
multiplication  and  division  will  produce  errors  in  the  fixed-point  arithmetic,  and  afi  elemen¬ 
tary  operations  will  introdue  errors  in  floating-point  arithmetic.  Therefore,  in  the  following 
we  discuss  the  two  case  seperately. 

For  fixed-point  arithmetic,  the  model  of  error  analysis  is  given  by 

fx(a  ±  b)  =  a  ±  6, 

fx(ab)  =  ab+  f , 

«!>-!  +  <• 

where  a  and  b  are  fixed-point  numbers.  £  and  (  are  errors  and  satisfy  |£|,  |£|  <  A,  where  A 
depends  on  the  configuration  of  the  computer.  Usually  A  =  C0~b,  where  /?  is  base  which 
computer  uses,  b  is  the  wordlength,  and  C  is  a  constant. 

For  floating-point  arithmetic,  the  model  is 

fl(a  +  6)  =  (a  +  6)(1  +  £)» 

fl(ab)  =  a  6(1  +  C), 

flfy  =  T{1  +  rt)' 

where  £,  t)  are  errors  and  also  satisfy  |£|,  |£|,  |t?|  <  A. 

Now  using  these  models,  we  study  the  effect  of  finite  wordlength  on  the  SLA.  We  denote 
the  computed  values  by  superscript  c  as 

3i,j  =  si,j  +  ai,j  (3.23) 
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4  =  4  +  /?,• 

tci  =  U+V 


Pi  =  + 


In  the  computer,  the  SLA  becomes 


(  \ 

/  \ 

4-i 

1 

4  = 

<  1  ) 

+ 

\  / 

4 


'  i  ' 


4-2 


\  1  / 


+  4 


4 


n 


t 


n-i 


~  +  T}i 

-i 

1  +  Qi- l4-l  +Ci  +  <Ti 


(3.24) 

(3.25) 

(3.26) 


(3.27) 

(3.28) 

(3.29) 


In  the  above  expression,  the  errors  0,,  77,,  <r,  are  local  errors  and  a.j ,  7,,  /x*  are  total 

errors.  Note  that  we  will  consider  only  the  first  order  error  terms  in  our  analysis;  the  second 
and  higher  order  terms  will  be  ignored. 

Now,  we  derive  the  relationship  between  local  and  total  error.  Using  (3.23),  (3.25)  and 
(3.29),  we  have 


7.-  =  n-u 

=  1  +  &-i4-i  +  Ci  +  <7i  -  (1  +  a-rli-i  +  Cj) 

=  Si-iSU-1  +  Vi  (3.30) 


Taking  the  same  procedure,  we  also  have 


*  jb  j  _  U  7 i  ~  47.-1  ,  _ 

Pi  -  di  ~  di  -  —  +  rji-  - —  - - - +  rji 

s_i  t.-i  if— 1 


(3.31) 
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where  £*i  =  sf  —  si  =  —  2t)i  and  a 2  =  —  2rji  —  772.  The  residual  is  defined  as 


(3.32) 


rn  —  Tnsn  Tnsn  —  Tnan 


(3.33) 


Substituting  (3.32)  into  (  3.33),  we  have 


(3.34) 


Since  T,  can  be  partitioned  as 


(3.35) 


(3.36) 
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Similarly,  partitioning  Ti  as 


we  get 


Also,  partitioning  Ti  as 


We  may  write, 


T,= 


'  .  d_, ' 

Qi-i 


Ti 


/  \ 
0 


(a  a.  ' 


Ti  = 


■  / 

\ 

1 

Qi— 2  C»— 1 

Oi-2 

Ti- 2  Qi-2 

Ci-l 

Cl  2  1 

<  0  ' 


Ti 


eu-2 

o 


) 

t  \ 

g-iQU-2 

Li-2 

\  4-2SU-2  ) 


Substituting  (3.36),  (3.38)  and  (3.40)  into  (3.34)  and  noting  that 

1  +  Si_2^.-2  +  °i- 1 

Oi-2  +  Ti-2§i-2  +  Ci-2 


t  \  ( 

1 


Ti 


Si-2 


\  1  ) 


*i-l 


\ 


Ci-l  +  C-_25j_  2  +  1  J 


Q 

V  u~x  ) 


(3.37) 


(3.38) 


(3.39) 


(3.40) 


(3.41) 


we  get 


(  \ 

t  \ 

Li- 1 

Cj-ifiU-i 

Li  = 

+ 

V  U~x  > 

(  0i-2SU-2  \ 
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+  Tiei 


(3-42) 


For  the  discussion  of  the  bound  of  the  residual,  we  seperate  (3.43)  into 


—<Tj  +  dj<Tj-i  —  TJjtj-i 

£-i,j  =  0  +  TjQ-j 

K  -ffi  +  djCTj-i  -  rjjtj- 1  j 

=  A 


(3.43) 

(3.44) 

(3.45) 

(3.46) 

(3.47) 
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Suppose  that  m  =  J2j=i  Ej,j  is  true  when  /  <  i  —  1.  Then 

(  v^«-l 


Efij  = 

i= l 


-  4 


£;=i  £.-x.i 
o 

/ 
l 

0 

o 


w  0  x 


+ 


e;=i  £i- 


/ 

Li- 1 

+ 

/  \ 
0 

<  0 

^  > 

( 

-a 

TiSd  + 

Q 

^  -(Ti  +  di<Tj_ 

-d, 


'  0  X 


£-2 


\  0  / 


=  Li 


(3.48) 


Note  that  r\  =  —2 (<tx  -f-  j/x)  +  Ti£a  =  F\t\.  Thus  (3.47)  is  true  by  mathematical  induction. 
Here,  the  vector  F_j,j  represents  the  error  produced  during  the  computation  at  jth  stage, 
and  Fj  j  represents  the  error  produced  at  jth  stage  propagated  into  s,-. 
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3.4  The  bound  of  residual 


We  have  obtained  the  recurrence  relation  of  the  residual  in  the  previous  section.  In  this 
section,  we  derive  a  bound  on  the  residual  vector.  We  know 

li-i  =  a.-i  +  o.-i  (3-49) 

Since  ||o1_1||  <  flmiiCl  +  |Pm|)  -  1  [6],  have 

llfc-ill  <  2  fn(l  +  |Pn.|)-ll  (3.50) 

.m=l 

We  will  use  this  relation  in  the  following  discussion.  Before  we  get  the  bound  of  residual, 
we  firstly  prove  two  lemmas. 


Lemma  5  If  a  symmetric  vector  u,-  satisfies  the  recurrence 


di  =  (1  +  p,)(l  -  pi- 1) 
then  another  vector  can  be  defined  as 

Ai  =  (1+  Pi) 

The  vector  Ai  satisfies  the  recurrence, 


(3.51) 

(3.52) 

(3.53) 

(3.54) 
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Ml  <  (2i  -  1)114,11 


(3.55) 


Proof:  Substituting  (3.53)  into  both  sides  of  (  3.54),  we  obtaine  (3.51).  This  prove  that  we 
can  define  A,  as  in  (3.53).  Noting  that 


Ai  =  (1  +  Pi) 


(  \ 

v*-i 


\  0  / 


—  Vj 


we  have 


A:  = 


(  A  \  (  -  \ 
'  -  pi 


\  0  I 


0 

V  4-i  j 


From  (3.53)  and  (3.57),  we  have 


(3.56) 


(3.57) 


vi  =  (Si-Ii)~1(Ai-Ai)-Ai 


(3.58) 


where 


/  \ 

0  0  •••  0  0 

1  0  •••  0  0 

0  1  •••  0  0 


(3.59) 


0  0 


1  0 


and  Ii  is  the  identity  matrix  of  dimension  i.  The  inverse  of  5,-  —  /,•  is  the  lower  triangular 
matrix  having  ail  elements  on  and  below  the  main  diagonal  equal  to  —1.  So  from  (  3.58), 
we  have 


&  =  (Si  -  Ii)~l£i  +  Li  Ad-, 


(3.60) 


where  L,  is  the  lower  triangular  matrix  having  all  elements  on  the  main  diagonal  equal  to 
0  and  all  elements  below  the  main  diagonal  equal  to  unity.  The  norm  of  is  bounded  as, 

fell  <  ||(A  -  All  dj AAll 
<  *1IAII  +  (*-  1)||A|| 

=  (2i  —  1)||  A  ||  (3.61) 


At  this  point,  we  have  proved  the  lemma.  Now  we  state  and  prove  Lemma  6. 
Lemma  6  If  a  vector  A,  satisfies 


4+1 


(  \ 

(  \ 

A 

0 

—  p*+i 

<  0  > 

<  A  > 

and  ||  A  |l  ^  B,  then 


iiAii  <  b  n  (i+i^D 

m=:j+l 


Proof:  Suppose  that  ||A||  <  5rim=i+i(l  +  l/»m|),  then 

IIA+ill  <  II A II  +  |p,+i|||A|| 
^  (!  +  l/i.+iDIIAII 
<  B  ff  (l  +  bml) 

HI  =  7  +  1 

Since  ||Aj||  <  B.  By  induction,  we  have  proved  the  lemma. 
Using  the  Lemmas  5  and  6  in  (3.46),  we  get, 


\\e.j\  <  (a  -  i)b  n  (i+w) 

m=j>'+l 


(3.62) 


(3.63) 


(3.64) 


(3.65) 
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where  B  is  determined  by  F_jj  as  defined  in  (3.44).  From  (3.44)  and  (3.53),  we  have 


4  = 


For  fixed-point  arithmetic,  we  have 


Wi\  <  j A 


3~  2 


<  2A  (1  +  | pT 

m=l 

l«ul  <  A 


and 


Il£j.jll  ^  2|<Tj|  +  2(f,|<Tj_l|  +  ||7>£j||  +  1\T),ty-l\ 

<  A(/  +  4)  nil  +  KI) 

m=l 

=  B 


Substituting  (3.70)  into  (3.65),  we  get 

ll£ij  <  (2i  -  1)A0'2  +  4)  ri  (1  +  M) 

m=l 

and 

lln II  <  EII&J 

j=i 

_  .  ,2i'  2i3  47i3  25i‘„  ,  .. 

-  A(T  +  T  +  1  6“  1  + 


(3.66) 


(3.67) 

(3.68) 
(3.60) 


(3.70) 


(3.71) 


(3.72) 


Writing  the  above  equation  for  i  =  n,  we  obtain  (3.21).  For  floating-point  axithmetic,  we 
have 


I  <  ;'A(2  +  ||5,_i||) 


(3.73) 


and  since 


so 


where  (,  £, 
as 


Mj-il 


<  2j a  n  (i + iP„)i 

m=l 

<  Alij.,1 

<  aafid  +  W) 

m=l 


(3.74) 


fKai)  = 


fi 


(  \  (  \ 

Si-i 

'  + 


L\  1  / 


l 

\  )\ 


-fl 


di 


/  \ 
1 


S.i-2 


\  1  /JJ 


=  (1  +  0 


(  \ 

(  \ 

' 

> 

(1+0 

Si-i 

+ 

1 

» 

K  1  > 

<  1‘-1  > 

4 

“  (1  +  0(1  +  v)di 


{  \ 
1 


Si-2 


\  1  / 


fli  =  (C  +  0 


/  \ 

*-l 

L\  1  / 


+ 


/  .  \ 


-  (C  +  n)di 


•  i  ' 


Si-2 

1 


(3.75) 


\  A  / 

7  axe  local  errors  and  satisfy  |£|,  |£|,  [77I  <  A.  The  bound  on  0 ,  can  be  written 


114,11  <  2A  [2  +  2||iJ_1||  +  2d,  +  d,||^_3||] 

<  12in(HW) 


(3.76) 


m=l 
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so 


ll£;J  <  2|ffjl  +  2<iiky-i|  +  llr^ll  +  2\Vjtj-i\ 

<  A(20i  +  4)n(l+M)-B 

m=l 

from  (3.65),  we  get 

ll£ij  <  (2i  -  l)A(20j  +  4)  n (1  +  M)  (3.77) 

m=l 

and 

IlCill  <  (20i3  +  18i2  -  14*) A  n  (1  +  \pm\)  (3.78) 

m=l 

Writing  the  above  equation  for  i  =  n,  we  obtain  (3.22). 

From  these  bounds  we  see  that  if  |/)j|  is  not  close  to  one  and  n  is  not  large,  the  residual 
is  small.  This  means  that  the  solution  of  Toeplitz  system  is  accurate.  If  n  is  large,  and  the 
problem  is  ill-conditioned,  the  residual  may  be  large. 

It  is  worthwhile  to  mention  here  that  we  cannot  define  A,  by  u,-  =  A,  +  A,  in  lemma  5  as 
this  definition  will  require  a  polynomial  devision  by  (1  —  z_1)  to  get  A*  from  un.  However, 
it  is  not  guaranteed  that  this  requirement  will  be  satisfied  due  to  errors.  In  filter  design, 
this  means  that  a  pole  exists  on  the  unit  circle  and  the  corresponding  filter  is  unstable. 
This  could  mean  that  the  SLA  is  not  strictly  stable.  Consequently,  we  have  to  adopt  an 
alternative  route  to  reach  our  objective. 
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Chapter  4 


The  Fast  and  Superfast  Algorithms 
for  Block  Toeplitz  Matries 


In  the  last  two  chapters,  we  have  studied  the  fast  algorithms  to  solving  Toeplitz  systems. 
The  algorithms  have  the  complexity  of  computation  with  order  0(n2).  They  are  more 
efficient  than  the  classical  methods.  However,  this  is  not  the  end  of  searching  for  more 
efficient  algorithms.  Another  category,  superfast  algorithm,  that  is  ,  algorithms  with  order 
0(ralog2n)  complexity,  will  be  studied  in  this  chapter. 

Until  now,  the  Toeplitz  system  that  has  been  studied  corresponds  to  the  scalar  case, 
that  is,  the  elements  of  Toeplitz  matrix  are  scalars.  The  block  Toeplitz  system,  that  is, 
each  element  of  the  Toeplitz  matrix  is  also  a  submatrix,  appears  in  applications  such  as 
image  processing.  The  fast  and  superfast  algorithm  for  the  block  case  is  also  studied  in  this 
chapter. 

We  firstly  study  the  fast  algorithm  'or  block  Toeplitz  system.  We  proved  the  Gohberg- 
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Semencul  formula  for  a  general  Toeplitz  matrix, and  then,  introduce  the  Schur  algorithm 
for  the  block  matrix.  The  superfast  algorithm  will  be  studied. 


4.1  The  Fast  Algorithm  for  Block  Toeplitz  Matrices 

The  Levinson  algorithm  for  solving  the  scalar  Toeplitz  systems  has  been  studied  in  the 
previous  chapters.  The  details  of  algorithms  for  inverting  a  block  Toeplitz  matrix  are 
studied  in  this  section  [16,  17,  18].  The  results  for  the  block  Toeplitz  matrix  are  very 
similar  to  those  for  scalar  case. 

Consider  a  Toeplitz  matrix 

Tp,n  =  (Ti-j)  i,j  =  1, 2,  •  •  • ,  n  (4.1) 


where  each  T,  is  p  x  p  matrix.  We  want  to  get  a  fast  algorithm  which  gives  an  inversion  of 
T 

p,n* 


We  can  partition  Tp,n  into 


where 


(4.2) 
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In  this  chapter,  superscript  t  represents  block  transpose,  and  ’  represents  the 
transpose. 

For  a  matrix 


We  know  the  inverse  of  A  is 


A'1  = 


An  A12 

A2i  A22 

on  0x2 
021  022 


where 


an  =  Axx1  +  Ajj1  Ax2(A22  —  A21A1X1  A12)  1A2iA{“11 
012  =  —  AJ"11Ai2(A22  —  A2iAj‘11Ai2)-1 
°21  =  —(A22  —  A2lA1i1Ai2)  1A2iA111 
«22  =  (A22  —  A2lA^x1A12)-1 

So  if  we  use  this  formula  for  Tp,n,  we  have 

T~l  ,  +W  a~lVl  W  a-1 

_1  p,n-l  ’  ii_nun  T—n  X!-nun 

^p,n  = 

OJ1 

where 

W  —  -T~l  T 

iJ-n  —  1  p,n-l— -n 

A  f  A  * 

K  =  -XJ'-U 

an  —  To  —  XnTp,n-lt.-n 


common 


(4.3) 

(4.4) 

(4.5) 

(4.6) 
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From  (  4.4  ),  we  get 


(4.7) 


(4.8) 


(4.9) 


-T'Lit,  (4.10) 

-t. (4.ii) 

To  -  (4.12) 

Substituting  (4.9)  into  (4.7),  we  have 

Kn+l  =  - 

where 

=  lnl-n  +  ^-(n+l)  (4.14) 
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In  this  chapter,  we  define 


Jp,n  — 


0 

0 


Note  that  Jp^n^Ip^n  —  Iptn* 

From  (4.7),  we  finally  have 


0  h 


Ip  0 


0  0 


J£n+1  = 


■  * 

•  * 

In 

2Ln 

— 

0 

.  ^  . 

A.  “t. 


(4.15) 


(4.16) 


From  (4.5),  we  have 


i-n+l  =  ~I-n+l 


T  T  T~* 

t'p,ntJp,nJ-  p,n— 1 


So 


where 


H£+i 


rpt  rp—t 
—JLn+l-lp<n-l 


=  TlJLn  +  T, 


n+1 


=  T* 

i-n+l 


Taking  the  same  procedure  on  2£n+l  and  2Ui.  we  have 


zUi 


•  * 

— n 

Kn 

0 

.  . 

• 

* 

r  —1 

*  t  . 

£ 

0 

—  fnan 

X*  Ip 

(4.17) 


(4.18) 


(4.19) 

(4.20) 
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where 


=  Zll»  +  Tn+1 

(4.21) 

=  zUiL + T_(n+1) 

(4.22) 

Note  that  the  un,  on,  en,  an  and  /3n  are  p  x  p  matrices.  Now  we  prove  some  of  relations 
between  un,  on,en,fn,an  and  /?„. 

Using  (4.4)  and  (4.11),  we  have 

«.-/.  =  Zll-n  -  E-jL  =  0 


or 


By  the  same  steps,  we  can  prove 


un  —  /n 


°n  —  e„ 


Let  us  find  the  relation  of  an  and  /3n  now.  From  (4.6),  we  have 


®M-1  -  ®n  =  ~On0~'u„ 


(4.23) 

(4.24) 


or 


similarly, 


From  (4.9),  we  have 


«n+l  =  <*n  ~  Onp~lUn 


/^n+1  /^n  fn&n 


T~l  =  J  T-t  7 

ip,»  JP,n1p,nJP,n 


(4.25) 

(4.26) 
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we  observe  that 


p-1  B-'ZL 

XA“x  r-^  +  Xn/5-1^ 


TO  oj  =  ft 


-1 

n 


h  zl 


(Tp,n)i+hj+l  ~ 


X, , 


TOi  +  JCA^ki 


From  (4.3),  we  also  have 


Cki-fti+Mlw 


(4.27) 


(4.28) 

(4.29) 

(4.30) 


(4.31) 


Combining  (4.30)  and  (4.31),  we  get 


(r,-,1  w»  =  (Oj  +  -  tiUci'&u 


The  equation  (4.32)  can  also  be  written  as 


T-i  _ 

P.n 


/p  0 

Xx  IP 


0 

0 


xn  •••  Xj  7p 


0  0  •  •  •  0 


ft*1  0 


0  A 


0 


W\  •••  Wn  0 


a 


-l 

n 


0 

0 


0  0 


-i 


-i 


0  0 


0  a"1  0 


0  0 


a, 


-l 


0  0 

0  V„ 

0  0 

0  0 


(4.32) 


Ip  Z\  •  •  •  Zn 
0  Ip  •  •  •  Zn_i 


Vi 

Vi 


(4.33) 
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This  is  a  very  important  result. 
From  (4.4)  and  (4.5),  we  have 


Now 


Similarly, 


Now  let 


Tl^Wn  +  Z-n  =  0 

(4.34) 

=  0 

(4.35) 

0 


(4.36) 


Ip  YL 

T*  = 

p,n 

an  0 

h  & 

T  — 

■Lp,n 

o 

Ip 

Pn 

2L 

0 

(4.37) 

(4.38) 

(4.39) 


£ 


U  £ 

Ip 

2£n 

ip  n 


\ 
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w„ 


From  previous  discussion,  we  have 


•*p 


TLnKn  = 


y‘  Tx  = 

— ™  J'Pin 


Tp,n2L  = 

&Tp,n  = 


W 


■n+l  — 


J£+i  = 


*n+l  = 


zUi  = 


Ctn 
0 

an  0 

fin 
0 

fin  0 

Kn 
0 

Zn  0 

X , 

0 

zi  o 


o 

In 

OnftT1 

0 

A 

i£„ 

-  /na;1 


fin 


A  f 

0  z 


«n_1«n 


0  £ 


or  in  polynomial  form 


(4.40) 

(4.41) 

(4.42) 

(4.43) 

(4.44) 

(4.45) 

(4.46) 

(4.47) 


Wn+l(2)  =  Wn(z)  -  zXn(z)0?Un 

(4.48) 

K*  (*)  =  VZ{z)-zonfi?Z'n{z) 

(4.49) 

Xn+i(.z)  =  Xn(2i)  -  zWn(z)a‘1en 

(4.50) 

=  %(*)  -  *fnd:iyi(z) 

'  (4-51) 
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where  WQ  =  Vo  =  X0  =  Zo  =  IP  and  aQ  =  @0  —  T0.  We  can  see  that  (4.48)-(4.51)  have  the 
same  recursive  form  as  the  Levinson  algorithm. 

If  we  set  the  p  =  1,  that  is,  Tp>n  is  a  scalar  matrix,  then  from  (4.6)  and  (4.12),  and 
noting  that  (3n  and  an  are  scalars,  we  get 

=  OCn 

From  (4.40)  —  (4.43),  we  have 

2L=Vn  and  Wn  =  Zn 

It  is  clear  that  these  are  two  unknown  vectors  in  this  case  and  the  algorithm  can  be  rewritten 
as 


Wn+i(z)  =  W„(z)  -  zu}n+iXn(z) 

(4.52) 

Xn+1(z)  =  Xn(z)  -  zin+xWn(z) 

(4.53) 

where 


w„+i  =  0~'un  =  a“Vn 
7n+l  =  On/?"1  =  allen 

This  is  the  generalized  Levinson  algorithm  for  a  scalar  Toeplitz  matrix  which  is  non- 
symmetric. 


If  we  assume  that  Tp<n  is  symmetric,  then  we  also  have 


and 


wn+ 1  =  7n+l 

The  algorithm  becomes  the  classical  Levinson  algorithm  which  was  studied  in  chapter  2. 


4.2  The  Schur  Algorithm  for  Block  Toeplitz  Matrices 


In  this  section,  we  use  the  transformation  on  the  result  of  last  section  to  get  Schur  algorithm 
on  block  case. 


Let 


Then 


En-k  Q 


-iLn—k 


Ok  In—k+1 


T 

±j>,n 


2Lk- i  0 
0  W'k-i 

o  o 


En-k  = 


Ok  In—k+1 


T 

±P,n 


Km-1 

0 


Ok  In—k+1 


I'ptk—l  * 

L\  * 


2Lk- 1 

0 

0 


=  LxKm-i 


(4.54) 


(4.55) 
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(4.56) 


The  first  row  of  L\  is  t!  so 


*  E 

Pn-kfi  =  Hk2Lk-l  =  °k- 1  =  eAr-l 


Note  the  star  *  represents  the  term  which  we  do  not  care  about.  For  Q_n_k  we  have 


0 


9-rx-k  - 


Ok  In-k+ 1 


L  Pt^ 


IKjb-i 

0 


Ojt  4i-fc+ 1 


Tp,*  * 

L2  * 


0 

±fc-l 

0 


SO  Qn-kfl  =  Ojt-1. 

Now  we  axe  going  to  obtain  the  recurrence  relation  of  and  Q_n_k- 


£Ln-k-l  — 


0k+ 1  In-k 


T 

■LP,n 


0 


(4.57) 
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Ok+ 1  In-k 


LP.n 


■ 

. 

> 

2U-i 

0 

0 

_ 

Kk-x 

: 

0 

k 

0 

. 

J 

Pn-k,  1 

Qn—k,  1 

■1° n—k,n—k 

Qn—k,n—k 

<*k-lek-l 


Another  recurrence  relation  of  Q  ,  is 

-2-n— k 


Qn-k-i  - 


0k+ 1  In—k 


LP,n 


0 

!£ 

0 


Ojt+i  /„_*■ 


TP,n< 


0 

0 

0 

2L- 1 

lL-i 

— 

0 

Pk-\Uk-l 

0 

1 

; 

0 

J 

Qn—k,  0 

P a—k,0 

Qn— fe,n— fc— 1 

Pn—k,n—k—l 

Pk-\uk-\ 


(4.58) 


(4.59) 
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The  initial  condition  is 


Q.n—1 


T,  To 


Tn  Tn-i 


Now  we  define 


Rn-k  Rn-k 


Ok  In~k+l 


■p,n 


Wk_  1  0 

o 

0  0 


Taking  the  same  procedure,  we  can  prove  that 


Rn-k-l  — 


Rn—k—l  ~ 


Rn-k,  1 


Rn—k,n—k 

Rn—k,  0 


Sn-kA 


Rn-k,n-h 


Pkhuk-i 


Sn—k,n—k—l 


Rn—k,  0 


Rn—k,n—k—l 


Rn-k,0  =  fk-1 

Sn-k,  0  =  #fc-l 


(4.60) 


(4.61) 


(4.62) 


(4.63) 


The  initial  condition  is 


Rn-l  Rn-\ 


r_i  To 


T_n  T_ 


n+1 
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If  we  write  the  equations  obtained  in  this  section  in  polynomial  form,  we  have 


zPn-k-\{z)  =  Pn-k(z)  -  Qn-k{z)ctZlxtk-\  (4.64) 

Qn-k-i{z)  =  -Pn-k{z)^k-iuk-\  +  Qn-k(z)  mod  zn~k  (4.65) 

zRn-k-i(z)  =  Rn-k(z)  -  Sn-k{z)^k\uk-l  (4.66) 

5n_fc_i(z)  =  -Rn-k(z)a.kl\ek-\+  Sn-k(z)  mod  zn~k  (4.67) 

where 

Pn-kfi  =  Cfc-lj  Qn-kfi  —  <*k-l  (4.68) 

■S'n-fc.O  —  ^k- 1>  Rn-k,0  =  «fc- 1  (4.69) 


4.3  The  Superfast  Algorithm  for  Block  Toeplitz  Ma- 

o 

trices 

In  the  previous  two  sections,  we  have  discussed  the  Levinson  and  Schur  algorithm  for 
block  Toeplitz  matrix.  Since  their  recursive  structure  is  same  as  Levinson  algorithm,  the 
complexity  of  these  algorithms  are  0(n2p3).  We  study  the  feist  Fourier  transform  (FFT) 
based  superfast  algorithm  in  this  section. 

Let  us  define 

flfc+i  =  0kluk 

r*+i  =  akxek 
A*+i  =  ok0ix 


(4.70) 

(4.71) 

(4.72) 
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Afc+i  =  /tat1 


(4.73) 


so  we  can  write  (4.48)-(4.51)  as 
Xk(z)  Wk(z) 


Xk-t(z)  Wk-i(z) 


1  -fife 
-zTk  z 


Z'k{z) 

1  -A  kz 

ZLi(z) 

.  . 

-Ajb  2 

and  (4.64)-(4.67)  as 

Pn_Jt_l(2)  Qn-k-i(z ) 

Rn-k-liz)  Sn-k-l{z) 

L 

Now  let 


Pn-k(z)  Qn~k{z) 


Rn-k{z)  Sn-k{z ) 


i  -n* 
}r*  i 

7  -r* 

i 


-in* 


Ak(z)  Bk{z) 

k 

-n 

i  -n, 

zCk{z)'  zDk(z) 

«=1 

— 2r,-  2 

so 


Xk(z)  Wk{z) 

We  also  can  write  (4.76)  and  (4.77)  as 

Pn-k-l(z)  Qn-k-l(z ) 


Rn-k-l(z)  Sn.k-l(z) 


Ip  Ip 


Ak(z)  Bk(z ) 
zCk(z)  zDk(z) 


Pn-l(z)  Qn-l(z) 


Rn-l(z)  Sn.i(z) 


k 

n 

i=i 


1  k 

n 

i=i 


_Ir. 

r1  • 


in, 


-n, 

i 

-r, 

i 


(4.74) 


(4.75) 


(4.76) 


(4.77) 


(4.78) 


(4.79) 


(4.80) 


(4.81) 
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Now  we  prove  the  relation 


z  kAk{z )  z^Bkiz) 

k 

=n 

\  -a 

_  z~kCk(z)  zl~kDk(z) 

t=i 

.  ~l*Ti  1 

(4.82) 


Proof:  For  k=l,  (4.82)  is  true.  Suppose  that  (4.82)  is  true  for  k  =  /,  then  for  k  =  /  +  1 


f+i 

i  -a 

1 

^  -a 

~  —Cli+i 

n 

2 

=  n 

2 

Z  ' 

»=i 

.  1 

1=1 

-}r,-  1 

-Kw  1 

z  1  Ai(z)  z1  lB[(z ) 
z~lC,{z)  z1~lDl(z) 


1  _a+i 

~7ri+i  1 


2:  1  (Ai(z)  -  zBi(z)rt+l)  z  1  (-Ai(z)nl+i  +  xB[(z)) 

z-1-1  (Ci(z)  -  zDt(z) r/+1)  (-Ctiz)^  +  zD,(z)) 


Using  (4.78),  we  can  prove 


/+i 

n 

1 

z 

— n, 

_ 

r-('+»/l,+1(Z)  z‘-<'+')S,+l(z) 

1=1 

.  -K« 

1 

z-('+>lC,+1(2)  z‘-('+1  >A+i(r) 

Thus  (4.82)  is  proved  by  induction. 

Using  the  same  procedure  we  can  prove 


Jk 

;  r< 

Wi)  c*(i) 

n 

1=1 

1  . 

Wi)  Mi) 

m  - 

Now,  we  let 


Ek(z)  zFk(z) 

k- 1 

-n 

1  -A  k-iz 

Gk{z)  zHk(z) 

1=0 

-A  k-i  z 

(4.83) 


(4.84) 


(4.85) 
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so 


(4.86) 


Zl(z)  Ek(z)  zFk(z) 

V2(*)J  [<?*(*)  zHk(z) 

At  this  point,  we  can  use  the  divide-and-conquer  technique  to  compute  An(z),Bn(z),  Cn(z), 
Dn(z),  En(z),  Fn(z),  Gn(z)  and  Hn(z). 

The  outline  of  the  superfast  algorithm  is  shown  below. 


1.  If  n  =  1,  then 


An{z) 

— 

1 

Bn(z) 

— 

Cn(z) 

= 

Qn-l,0^n-l,0 

Dn{z) 

= 

1 

En{z) 

— 

1 

Pn(z) 

= 

~Rn-l,oQn-l,0 

Gn(z) 

= 

~Pn- l,oSn-l,0 

Hn(z) 

= 

1 

2.  If  n  >  1,  then  call  the  algorithm  to  compute  Az(z),Bz(z),  Ca(z),  I?a(z),  E%(z), 
E$(z),  G%(z)  and  Hn{z). 


3.  Compute  P£_i(z),Qf-i(z),  R^-i(z),S^.-i(z).  From  (4.80)-(4.82)  and  (4.84),  we  have 


Pf-i(z)  Q$-i(z) 


Pn-l(z)  Qn-l(z) 


z~%Aa(z)  zx~*B*(z) 
z~$Cn(z)  zl~?D^{z) 
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Rn-l(z)  Sn-i(z) 


tb%{\)  a%(\) 


4.  Call  the  algorithm  to  compute  A\(z),Bn. (z),  CS(^),  Z?a(z),  (2),  Fn(z),  Gn(z)  and 


^2(2). 


5.  Compute 


An(z)  Bn(z) 

A%(z)  B*(z) 

A%{z)  BUz) 

zCn{z)  zDn(z) 

zCf(z)  zDz(z) 

zCn(z)  zDn(z) 

En(z)  zFn(z) 

EUz)  zFn(z) 

2  2  v 

Ez(z)  zF*(z) 

Gn(z)  zHn(z) 

G*(z)  zHn.(z) 

2  2  '  '  J 

G$(z)  zH^(z) 

6.  Compute 

Xn(Z)  =  i4B(2)  +  2C„(2) 
Wn(z)  ~  Bn(z)  +  zDn(z) 
Z<(z)  =  En(z)  +  zFn(z ) 
VM  =  Gn(z)  +  zHn(z) 


Using  (4.33),  we  can  get  the  inversion  of  block  Toeplitz  matrix.  In  the  algorithm,  we 
use  the  FFT  algorithm  to  compute  the  multiplication  of  two  polynomials. 

Studying  the  algorithm,  we  find  that  each  stage  of  the  algorithm  requires  32  fast  Fourier 
transforms  (  we  take  the  inverse  FFT  same  as  the  FFT  ).  Since  each  block  is  p  x  p  ,  so  we 
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need  a  total  of  32p2  FFT.  The  number  of  arithmetic  operations  is 

M(n )  =  2A/(^)  +  32p2^log2n  +  0{p3n) 

=  &p2nlogln  +  0(p2nlog2n)  +  0{p3n) 

A(n)  =  2A(^~)  +  32p2nlog2n  +  0(p3n) 

=  \§p2nlog\n  +  0(p2nlog2n)  +  0(p3n). 

Whre  M{n )  is  the  number  of  MULT  and  A(n)  is  the  number  of  ADD  required  by  the 
algorithm. 

There  are  two  case  which  will  simplify  the  computation. 

1.  The  matrix  Tp,„  is  symmetric,  i.e.  Tp<n  =  Tp  n.  In  this  case,  from  (4.6)  and  (4.12),  we 
find  that 


Taking  transpose  of  (4.41)  and  (4.43),  we  have 


KnOC.)'  - 


Of 


and 


TP,  n(&)'  = 


0n 

0 


Comparing  them  with  (4.40)  and  (4.42),  we  get 


(&)'  =  UCi 
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and 


(£)'  =  *n 

2.  The  matrix  TPi„  is  symmetric  and  each  submatrix  is  Toeplitz,  i.e.  TPttl  is  Toeplitz- 
block-Toeplitz  matrix. 

This  case  has  been  studied  in  [18].  The  result  is 

On 

.  TP,n2L  = 

0 

A  A 

We  can  prove  in  this  case,  Ak(x)  =  Dk(x )  and  Bk(x)  =  Ck{x).  The  total  number  of 
arithmetic  operations  is 

M(n)  =  2.bp2nlog\n  +  0(p2nlog2Ti )  +  0(p3n) 

A(n)  =  bp2nlog\n  4-  0(p2nlog2Ti)  +  0(p3n ) 


Note  in  [18],  Jp,n  is  defined  as 

0  0  •••  0  Jp 

0  0  •••  Jp  0 

Jp,n  — 

Jp  0  •••  0  0 

If  Tp,n  is  not  symmetric  but  Toeplitz-block-Toeplitz,  we  have 


Tp,n2Ln 
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The  total  operations  are 


M(n )  =  k.bp2nlog\n  +  0(p2nlog2n )  +  0(p3n ) 

A(n)  =  9p2nlog\n  +  O(p2nlog2n)  +  0(p3n) 

4.4  Conclusion 

In  this  chapter,  we  studied  the  superfast  and  fast  algorithms  for  computing  the  inverse 
of  block  Toeplitz  matrices.  The  results  are  valid  for  any  kind  of  block  Toeplitz  matrix. 
Some  simplifications  axe  also  studied.  The  results  are  very  efficient  when  p  -C  n.  It  should 
be  mentioned  that  Toeplitz-plus-Hankel  matrix  can  be  transformed  into  Toeplitz-block- 
Toeplitz  matrix  [19].  In  this  case,  p  =  2.  So  the  superfast  algorithm  will  be  very  efficient 
in  this  case.  For  the  scalar  Toeplitz  matrix,  de  Hoog  has  studied  the  superfast  algorithm 
[10].  In  our  case,  p  =  1  corresponds  to  the  scalar  case. 
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Chapter  5 


The  Finite  Precision  Analysis  of 


Classical  and  Split  Algorithms 


The  Levinson  algorithm  is  widely  used  in  the  linear  prediction  because  of  its  simple  recur¬ 
sive  and  efficient  structure.  There  are  other  algorithms  that  have  also  been  used.  These 
algorithms  are  Schur  algorithm  and  lattice  algorithm.  Recently,  the  redundancy  of  these 
algorithms  have  been  shown  in  [9].  In  [9],  a  set  of  new  algorithm  defined  as  split  Levin¬ 
son  algorithm,  split  Schur  algorithm  and  split  lattice  algorithm  have  been  proposed.  The 
difference  between  classical  and  their  split  correspondent  is  that  the  split  algorithms  save 
computation  time  by  almost  50%. 

In  spite  of  the  wide  use  of  those  algorithms,  a  very  important  aspect  has  not  been 
studied.  This  aspect  is  the  effect  of  finite  precision  of  computer.  This  problem  is  inevitable 
in  the  implementation  oe  the  algorithms  because  of  the  finite  word-length  of  computer. 

The  finite  precision  of  algorithms  is  very  important  in  the  practical  situations.  For 


69 


example,  a  designed  filter  that  is  stable  in  theory,  may  becomw  unstable  after  computations 
by  the  computer. 

So  fax,  some  results  about  Levinson  and  lattice  algorithm  were  reported.  In  [20]  and  [14], 
the  authors  showed  that  implementing  Levinson  algorithm  using  finite  precision  arithmetic 
results  in  instability.  G.  Cybenko  analyzed  the  finite  precision  effect  on  the  stability  of 
Levinson  algorithm.  But  he  did  not  give  a  formula  to  compute  the  errors  of  reflection  coef¬ 
ficients  [6].  The  recently  report  on  finite  precision  effects  on  Levinson  and  lattice  algorithm 
was  done  in  [12].  They  give  two  simple  formula  for  computing  errors  due  to  quantization. 
However  ,  they  did  not  analyzed  other  algorithms  mentioned. 

In  this  chapter,  we  study  the  finite  precision  effects  on  Levinson,  Schur,  lattice  and  their 
split  cases.  Some  simple  but  useful  formulas  will  be  given  for  estimation  the  errors  in  the 
refletion  coefficients. 


5.1  The  Classical  Algorithms 

5.1.1  The  Levinson  Algorithm 

In  the  previous  chapter,  we  have  studied  Levinson  algorithm.  For  convenience  of  the  reader, 
we  will  briefly  reintroduce  it.  For  a  sequence  of  real  signal  samples  u(0),  u(l),  •  •  •,  u(n  —  1) 
of  finite  length  n,  the  autocorrelation  axe  defined  as 

°i  =  £  “(*)«(*  -  0 

t=o 
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The  k-th  reflection  coefficient  can  be  computed  by  using  the  Levinson  algorithm.  The 
procedure  of  computation  is  as  follows. 

For  k  =  1, 2,  •  •  • ,  n 


k- 1 


9k  =  Ck-iOk-l,i 

(5.1) 

:'=0 

9k 

Pk  = 

(5.2) 

<*k-i 

Ok,i  =  +  pkdk-l,k-i 

(5.3) 

Clk  =  Ctk-1  +  pk9k 

(5.4) 

The  initial  conditions  are  ao,o  =  1,  and  ao  =  Co.  In  these  formula,  is  a  predictor 
coefficient,  and  a*  is  a  prediction  error. 


5.1.2  Schur  Algorithm 


Schur  algorithm  can  be  derived  from  the  Levinson  algorithm.  This  has  been  studied  in  the 
last  chapter.  In  the  following,  we  will  reintroduce  the  scalar  Schur  algorithm  which  will  be 
analysed  in  the  later  portions  of  the  chapter. 

In  the  last  chapter,  we  obtained  the  general  Schur  algorithm,  which  is  valid  for  a  block 
general  Toeplitz  matrix.  Since  in  this  chapter,  we  are  only  dealing  with  the  symmetric 
scalar  Toeplitz  matrix  case,  the  general  results  can  be  simplified. 

Defining  the  vectors  pn-k,  9n-k  by 


Pn—k  Qn—k 


Ok  o 


Ok 


o  a* 


(5.5) 
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and  computing  £n_fc  and  ^n_k,  we  have 


k 


Pn—k—l,j  =  ^2ck+i+j-iak,i 

(5.6) 

1=0 

k 

9n— fc— l,i  =  ci+iak,i 

(5.7) 

i=0 


Introducing  the  polynomials 

n—k 

Pn-k(z)  =  ^2  Pn-k,iZ' 

t=0 

n—k 

qn-k(z)  =  ^2  9n-k,iZ* 

i=0 

and  taking  the  same  procedures  as  those  in  last  chapter,  we  obtain  the  recurrence  relation 

zpn-k-l(z)  =  Pkqn-k(z)  +pn-k(z)  (5.8) 

qn-k-i{z)  =  qn-k(z)  +  pkpn-k{z)  (5.9) 

Pn—k,  0 

Pk  =  - 

9»-fc,  o 

In  this  chapter,  we  will  introduce  one  set  of  variables  instead  of  pjn_k  and  <£n_k  by  defining 

k 

hk,j  =  5Zc|j-«|afe,»  (5<H) 

i=0 

where  k  +  1  —  n  <  j  <  n. 

Comparing  with  (5.6)  and  (5.7)  ,  one  can  note  that 


(5.10) 


hk,-j  ~  <ln-k-l,j 
hk,j+k+l  —  Pn—k—l,j 


The  Schur  algorithm  with  hk,j  is 

Pk 


hk-i,k 

hk- 1,0 


(5.12) 
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hk-i-j  +  pkhk-i,k+j 


(5.13) 


hk-j  — 

hk,k+j+ 1  =  Pkhk-\,-j-\  +  hk-i,k+j+l  (5-14) 

5.1.3  Lattice  Algorithm 

Lattice  algorithm  [21]  is  different  from  the  first  two  algorithms.  In  Levinson  and  Schur 
algorithms,  we  first  calculate  the  autocorrelation  coefficient  of  the  data,  then  solve  the 
Toeplitz  matrix.  In  lattice  algorithm,  we  directly  compute  the  reflection  coefficients  and 
forward-backward  prediction  errors  by  using  the  signal  samples. 

The  lattice  algorithm  can  be  implemented  by  a  filter  called  lattice  filter.  This  filter  is 
widely  used  in  the  speech  prediction.  It  is  also  the  basic  frame  for  ARM  A  model  [11]. 

The  two  important  terms  in  the  lattice  algorithm  are  the  forward  prediction  error  / 
and  the  bachward  prediction  error  6 k.  For  the  given  signal  u0,  m,  •  •  *,  u„_i,  we  can  form 
the  forward  prediction  error  of  order  k  or  bachward  prediction  error  of  order  k.  and  the 
forward  prediction  error  is  defined  as 

k 

fk,j  =  53  ak,iUj-i 
i~0 

The  backward  prediction  error  is  defined  as 

k 

bk,j  =  53  ak,k-iuj-i 

is  0 

Multiplying  by  on  both  side  of  Levinson  algorithm  in  (5.3)  and  the  reciprocal 
Levinson  algorithm  and  summing  the  result  from  i  =  0  to  i  =  k,  one  gets  the  lattice 
algorithm. 
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The  lattice  algorithm  is  summarized  as  follows. 


For  k  =  1,2,  ■  •  •  n 


n+k-2 


9k  —  )  ^  fk—l,ibk—l,i—l 

(5.15) 

i=l 

9k 

Pk  = 

<Xk-l 

(5.16) 

fk,i  =  fk-l,i  +  pkbk-l,i-l 

(5.17) 

bk,i  =  Pkfk—\,i  +  bk-ij-i 

(5.18) 

The  initial  conditions  are  /ot,  =  6o)t-  =  u,  and  cto  =  ul- 

We  also  can  express  the  lattice  algorithm  in  the  polynomial  form.  Let  us  introduce  two 
polynomials  which  are  defined  as 

n+Jfc-l 

/n+fc-l(*)  =  fk,iZ{ 

«=0 

n-ffc—l 

bn+k— ~  )  1  bkjZ 

i=0 

Substituting  the  polynomials  into  lattice  algorithm  (5.17)  and  (5.18),  we  have 


fn+k-l(z )  =  fn+k-2(z)  +  Pkzbn+k-2(z)  (5.19) 

bn+k-l(z)  =  pkfn+k-2(z)  +  zbn+k-2{z)  (5.20) 


The  initial  polynomials  axe 

bn-i(z)  =  fn-i(z)  =  uizn~1 

i=0 
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5.2  Split  Algorithms 


5.2.1  Split  Levinson  Algorithm 

In  the  previous  studies,  we  introduced  the  Levinson,  Schur  and  lattice  algorithms.  They  are 
called  the  classical  algorithms.  The  computational  complexity  of  these  classical  algorithms 
is  0(n2). 

Recently,  more  efficient  algorithms  were  proposed  [9].  These  are  named  split  algorithms. 
The  main  idea  of  split  algorithms  is  to  reduce  the  redundancy  existing  in  the  classical 
algorithms.  This  can  be  done  by  introducing  a  set  of  symmetric  or  antisymmetric  sequence. 
As  we  know,  every  function  is  a  combination  of  its  symmetric  and  antisymmetric  parts. 
Attributed  to  the  beautiful  structure  of  Toeplitz  matrix,  the  Toeplitz  problem  can  be  solved 
only  by  the  symmetric  or  antisymmetric  parts.  Since  the  only  half  the  coefficients  need  to 
specify  the  problem,  the  complexity  of  split  algorithms  will  be  reduced  by  one-half  time 
compared  to  their  corresponding  classical  algorithms. 

we  have  met  the  split  Levinson  algorithm  in  Chapter  3.  In  the  following,  from  the 
polynomial,  we  will  study  some  details  of  the  algorithm. 

Let  us  define  the  symmetric  polynomial  by 

»k(z)  =  ajt-i  +  z~xhk-i{z)  (5.21) 

The  polynomial  sk(z)  are  symmetric  since 

h  =  z~k[ak.1(z~1)  +  zak-i{z~l)\ 
r  z~xak-\{z)  +  afc_!(z) 
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Note  that  5^(2)  is  a  special  ak(z)  which  corresponds  to  the  reflection  coefficient  pk  equals 
one.  Under  this  situation,  the  Toeplitz  matrix  is  singular. 


The  antisymmetric  polynomial  tk(z )  can  be  defined  as 

tk(z)  =  ak-i(z)  -  z~xak.x(z) 

tk(z)  is  antisymmetric  since 

t(z)  =  z~k[a,k-\(z~l)  —  za.k-\(z~x)} 

=  z~lak-x{z)  -  ak.x{z) 

=  ~tk(z) 

We  will  focus  our  attention  to  the  symmetric  case. 

The  Levinson  algorithm  can  be  written  in  polynomial  form  as 

ak{z)  =  ak+i(z)  +  z-'pkak-liz) 

The  reciprocal  polynomial  of  ak(z)  derived  from  (5.22)  is 

ak(z)  =  z~lak-i(z)  +  pkak-i{z) 

Substituting  for  ak(z)  from  (5.22)  and  0^(2)  from  above  in  (5.21),  one  obtains 

(1  4-  pk)sk(z)  =  ak(z)  +  ak(z) 

Incrementing  the  index  in  (5.21),  one  gets, 

Sfe+!  (2)  =  ak(z)  +  z~lak(z) 


(5.22) 


(5.23) 
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Multiplying  both  side  of  (5.23)  by  z~l  and  taking  the  difference  with  s*:+i(z),  we  get, 

(1  -  z~1)ak(z)  =  Sjfc+x(z)  -  (1  +Pk)z~1sk(z)  (5.24) 

At  this  point  we  heve  established  the  relation  between  the  a.t(z)  and  s*;(z).  By  using  those 
relations,  we  can  derive  a  recurrence  on  s*;(z)  instead  of  the  Levinson  algorithm. 

The  derivation  of  the  recurrence  of  Sk(z)  as  same  as  the  derivation  of  stepdown  procedure 
in  Chapter  2.  We  directly  write  the  relation  as 

sk+i(z)  =  (1  +  z~1)sk(z)  +  /3fcZ-1Sjfc_i(z)  (5.25) 

where  /3k  =  (1  —  pk)(  1  +  Pk-i)  The  initial  polynomial  of  Sk(z )  are 

s0(z)  =  2,  Si(z)  =  1  +  z-1 

The  split  Levinson  algorithm  also  can  be  written  as 


k 


Tk  = 

(5.26) 

i=0 

A  =  — 

(5.27) 

Tk- 1 

■Sfc+l.i  =  sk,i  +  Sk,i-1  ~  /3kSk-l,i-l 

(5.28) 

.  A 

Pk  ~  1  .  , 

1  +  Pk-\ 

(5.29) 

with  the  conditions  Sk,o  =  1,  {k  >  1),  So.o  =  2. 

5.2.2  Split  Schur  Algorithm 

The  split  Schur  algorithm  can  be  dericed  from  the  Schur  algorithm  by  defining 

vkj  =  hk-i-j  +  hk-  i,k+i  (5-30) 
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or 


Vn-k(z)  =  qn-k(z)  +  Pn-k(z)  (5.31) 

where 

Vn-k(z)  =  Vk,0  +  Vk,lZ  H - +  Vk,n-kXn~k 

From  (5.31),  we  have 

(1  +  pk)yn-k(z)  =  Qn-k(z)  +  pkPn-k{z)  +  pkqn-k(z )  +  pn-k(z) 

Using  Schur  algorithm  in  (5.8)  and  (5.9),  we  have 


(1  +  pk)yn-k{z)  =  qn-k-i(z)  +  zpn—k—i (z)  mod  2rn  k  (5.32) 

Reducing  the  index  k  in  (5.32)  by  one,  one  have 

(1  +  pk-i)yn-k-i(z)  =  qn-k(z )  +  zpn.k{z)  mod  xn~k+1  (5.33) 

Substracting  (5.31)  by  (5.33),  and  substracting  zyn„k(z)  by  (5.33),  we  have 

yn-k(z)  -  (1  +  pk-i)yn-k+i{z)  =  (1  -  z)pn-k{z)  mod  zn~k+1  (5.34) 

zyn-k(z)  -  (1  +  pk-i)yn-k+i{z)  =  (z  -  l)qn.k(z)  mod  zn~k+1  (5.35) 


At  this  point,  we  are  going  to  derive  the  split  Schur  algorithm  from  Schur  algorithm. 
We  have  known  Schur  algorithm  by  (5.8)  and  (5.9).  Multiplying  (1  —  z)  to  (5.8),  we  have 

(1  -  z)qn-k-i{z)  =  (1  -  z)qn-k(z)  +  (1  -  z)pkpn-k{z )  mod  zn_fc+1 

Using  (5.34)  and  (5.35)  in  the  results,  we  have 

-zyn-k-i(z)  +  (1  +  pk)Vn-k(z)  =  -zyn-k(z)  +  (1  +  pk-i)yn-k+i{z) 

+Pk{yn-k(z)  -  (1+  pk-i)yn-k+\(z)  mod  zn~k+1 
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rearranging  the  terms  on  both  side,  we  get 


-zyn-k-i(z)  =  -(1  +  z)yn.k(z)  +  (1  -  pk)(  1  +  pk-i)yn-k+i(z)  mod  zn  k+1 

Finally  we  can  obtain  the  Schur  algorithm  after  removing  the  factor  z  on  both  side  of 
the  last  equation.  The  Schur  algorithm  has  the  form 


yn-k-i(z)  =  (1  +  z  )yn-k{z)  - /3kz  1yn-k+i(z)  mod  z 


n—k 


(5.36) 


where 


0k  =  (1  ~  />ifc)(l  +  Pk-\) 

-  Vk’° 

Vk- 1,0 

The  initials  of  the  polynomial  are 


Vn(z)  =  c0  +  y£22ciz  * 

»=i 

n— 1 

Vn-\{z)  -  XXc,-  +  C,+1)2r_1 

i=0 

We  also  write  the  Schur  algorithm  in  the  form  as 
For  k  =  1, 2,  •  •  •  n 


A 

Pk 


Vk,o 

Vfc-1,0 

1  - 


A 


1  +  Pk- 1 


=  VkJ  +  Vk,j+ 1  -  Awfc-lj+l 


(5.37) 

(5.38) 


(5.39) 

(5.40) 


vk+lJ 


(5.41) 


5.2.3  Split  Lattice  Algorithm 


Taking  the  same  procedure  as  in  previous  sections,  we  can  derive  the  split  lattice  algorithm 
from  the  lattice  algorithm. 

Let 

u>k,j  =  fk-ij  +  (5.42) 

or 

wn+k-i(z)  =  fn+k-2(z)  +  zbn+k-2(z)  (5.43) 

Using  the  lattice  algorithm  in  (5.19)  and  (5.20),  multiplying  (1  +  pk)wn+k-i(z),  we  have 

(1  +  p^Wn+idiz)  -  /„+*_! (z)  +  (z)  (5.44) 

Incrementing  the  index  k  in  (5.43)  by  one,  and  subtracting  it  with  (5.44),  we  have 

wn+k(z)  -  (1  +  Pk)wn+k~i(z)  =  (z  -  l)6n+fc_i(z)  (5.45) 

Subtract  (5.43)  with  incrementing  k  by  one  from  (5.44),  we  have 

wn+k(z)  -  (1  +  Pkzwn+k-i(z)  =  (1  -  z)fn+k-i(z)  (5.46) 

These  relations  built  up  the  bridge  between  wn+k(z)  and  bn+k(z)  and  fn+k(z).  They  will 
be  used  to  establish  the  split  lattice  algorithm. 

The  split  lattice  algorithm  can  be  obtained  by  multiplying  (5.19)  by  (1  —  z)  and  using 
(5.45)  and  (5.46).  The  result  is  shown  below 

Wn+k{z)  ~  (1+  pk)zWn+k-\{z)  =  Wn+k-\(z)  ~  (1  +  Pk-\)zWn+k-2(z) 

+PkZ((l  +pk-i  )wn+k-2(z)  -  wn+k-i{z)) 
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Simplifying,  it  we  have  the  split  lattice  algorithm  as 


wn+k(z)  =  u>„+fc-i(z)(l  +  z)  -  Pkzwn+k-i(z)  (5-47) 

where  =  (1  —  pjt)(l  +  Pk-i)-  The  initial  polynomial  are 

n— 1 

Wn-l(z)  = 

»=0 

n 

Wn(z)  =  ^(«i  + 

i=0 

The  split  lattice  algorithm  can  be  also  written  in  the  form  as 
For  k  =  1, 2,  •  •  ■  n, 


n+k— 1 

2t*  =  J2  w\  - 

(5.48) 

t=0 

A  =  — 

(5.49) 

^jt-i 

i  A 

*  =  1  i+*-, 

(5.50) 

wk+u  =  Wk,i  +  Wk,i- 1  -  PkWk-l,i-l 

(5.51) 

The  initial  conditions  are 

w0,i  =  2  u, 
lUl,.  =  Ui  +  Ui-l 


5.3  Finite  Precision  Analysis  of  Classical  Algorithms 

5.3.1  Finite  Analysis  of  Levinson  Algorithm 

The  finite  wordlength  of  computer  gives  rise  to  the  need  to  analyse  the  precision  of  output 
by  computer.  Due  to  the  finite  wordlength  of  computer,  some  errors  are  caused  in  the 
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results  obtained  by  performing  numerical  computations.  If  the  errors  grow  too  large,  the 
algorithm  can  not  be  used. 

For  analysing  the  finite  precision  effects  of  computer  on  the  results,  we  should  explore 
the  structure  of  algorithms.  Different  algorithms  will  react  in  a  different  manner.  Some  of 
them  give  up  a  precise  enough  result;  some  of  them  not.  For  the  analysis,  we  first  establish 
the  model  used  in  the  chapter.  Let  x  represent  a  precise  variable,  then  after  quantization 
it  will  become 

x  =  x  +  ex 

where  the  x  represents  the  approximation  of  the  variable  x,  and  ex  is  the  quantization  error. 

In  this  report,  we  only  study  the  fixed  point  arithmetic.  In  this  case,  since  the  output 
of  two  b-bit  number  multiplying  or  dividing  is  a  2b-bits  number,  the  quantization  errors 
will  occur.  For  the  addition  and  subtraction,  there  is  no  error  created  under  the  condition 
of  no  overflow. 

After  quantization,  the  Levinson  algorithm  becomes 


k- 1 

9k  ~  }  '  Cfc—jOfc—i,,’  *f-  £g,k 

i=0 

(5.52) 

9k 

Pk  =  -  ep,fc 

Olk-l 

(5.53) 

Qkji  =  •  lti  *4"  P  k^  k—\  tk—i  "1“ 

(5.54) 

d*  =  1  +  pkgk  +  Ca,k 

(5.55) 

with  the  conditions  a0,o  =  1  and  a0  =  cq. 

The  errors  of  the  computation  axe  obtained  by  subtraction  (5.1)-(5.4)  from  (5.52)-(5.55). 


For  example,  taking  (5.1)  away  from  (5.52),  we  get 


A  gk  =  gk  —  9k 

k- 1  fc-1 

=  51  Ck-i&k-U  +  eg>k  —  ^2  Ck-iak-l,i 

«= 0  1=0 

fc-1 

=  y  1  (cfc-.ajfe-i^  Ck—iQk—l,i  “t"  Ck—iQk—l,i  ck-iak-l,i)  4"  e3,fc 

t=0 

fc-1 

=  51  c^-«  +  Cfc-«A  flfc-i.t)  4-  e3ifc 

t=0 

fc— 1 

=  51  («fc-i,.A  cfc-«  +  Cfc-iA  afc-1,,)  4-  e3,fc  (5.56) 

«=0 

Where  we  have  defined  A()  being  an  operator  as  A  x  =  x  —  x. 

In  the  last  step,  we  use  ajt_iitA  c*_,-  to  replace  a^.^A  Because  ak-\,i  =  <ik-i,i  + 
Aan,,-,  so  a.k-i,iAck-i  =  ajt-i.i A  c*-,-  4-  A  afc_i,;A  The  last  term  are  the  product  of 
two  operator  A,  that  is  ,  A  ()A  ().  We  call  such  terms  second  order  errors.  In  this  report, 
we  assume  that  the  second  or  higher  order  errors  are  much  smaller  than  the  first  order,  and 
therefore,  they  are  ignored. 

Taking  similar  procedure  on  the  other  quantities  in  Levinson  algorithm,  we  have 


A  pk 


h  ~  Pk 

pkA  afc-i  +  A  gk 
afc-i  4-  Aafc_i 


Since  <**_!  Aaj._i,  we  can  simplify  the  result  as 


A  pit  = 


PkA  Qfc-i  +  A  gk 

Ok- 1 


eP,k 


(5.57) 


A  a*t,  =  a k,i  -  ak,i 

=  A  dk-i,i  4-  pfcA  afc_i,fc_j  H  ak-i,k-i  A  pk  4-  ea<k,i  (5.58) 
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Repeating  the  relation  Aa^  to  Aajk-i,,-,  Aajfc_2,t,  •  •  •,  we  have 


k 

A  ak'i  =  ^  (pj  A  flj-ij-t  +  ai-i,i-iA  pj  +  ea,j,i) 

i= i 

For  a*,  we  have 

A  a*  =  dk-  Ok 

=  Actk-1  +  PkAgk  +  9kA  pk  +  ea<k 

From  (5.57),  we  have  known 

A  £*  =  — afc_iApjt  —  pfcAa*,_i  —  afc_iePifc 
Substituting  A into  (5.60),  we  have 

A  afc  =  (1  —  p2)  A  arjt-i  —  PkOtk-x  A  />*  -  pkOCk-iePtk  +  gk  A  /s>jt  +  ea,fc 
Noting  that  —pkOk-i  =  <7/t,  we  have 

A  afc  =  (1  -  ^fc)A  a*_j  +  2yjt A  pfc  +  gk£p,k  +  eQiJt 
Continuing  to  expand  A  ajt-2,  A  a*_3,  •  •  •,  A  ao,  we  have 

fc  k—l  fe-l 

A  Ctk  =  JJ(1  —  Pi  0:0  d"  ^  n(l-  [2<jft'A  +  Cot.i  +  9i£p,i] 

i=l  i=0  j=i 

+2gk  A  />*  +  ea,fc  +  gkePlk 

Since  aj  =  n?=i(l  —  pl)ao »  we  can  simplify  Act*  into 

A  a*  =  — A  a0  +  —  [2^,-  A  pi  +  ea,»  +  9\*p,i[  +  2^*  A  pk  +  ea<k  +  gk^p,k 

ao  a* 

—  —  A  a0  +  oik  ^2  (-  2pi  A  pi  -  pit„,i  + 
ao  i^Q  \  ai  / 


(5.59) 


(5.60) 


(5.61) 


(5.62) 
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Substituting  A  a*  into  A  pk  as  in  (5.57),  we  finally  get  the  error  A  pk  as 

&Pk  =  -pk - -  +  Pk  ^2  +  PieP,i  -  — )  -  -  ep,k  (5.63) 

&0  i= o  v  OCi  J  OCk- 1 

At  this  point  we  disucss  the  result  of  (5.63).  The  complex  result  can  be  simplified.  The 
simplest  result  on  the  error  can  be  obtained  if  we  neglect  the  term  A  <**_!,  es>jt,  ePifc,  e0>jit- 
in  all  derivation  and  Acij_ij_t-  in  (5.59),  we  find  out  that 

Agk 


A  pk  = 


<*k- i 
1 


fc-i 


)  ]  (  ak— l,i A  Ck—i  ■+■  Ck—i  'y  Clj  —  iJ—i&Pj 

Olt-l  «=o  V  i= 1  , 


(5.64) 


This  formula  is  same  as  the  result  given  in  [12]. 

For  more  accuracy  we  can  continue  to  expand  A  ajt,,  as 

k  j- 1 

'A  Qfc,«  =  ^2  P'  ^2  (Pi\^ai\-i>i\-i+i  "b  aii -i.ii -i+i^Pii  "b  eo.ii.i-«) 
i- 1  ii=i 
k 

-b  ^2  Pj  +  ea,j,t) 

J=1 

Continuing  expansion  of  Aa*^,  we  find  that  A  a*,,,  is  a  first  order  function  of  terms  , 
pjPjx ,  •  •  •  and  A  pj ,  that  is 


where 


A  =  5Z  /(/>;>  PjPh  >  •  •  •) A  ft 

i=x 


/(ft)  Pi  Pin’  ”)  —  fo  +  fiPj  +  flPiPix  + 


(5.65) 


If  we  ignore  the  terms  pj,  PjPjx ,  •  •  •  and  A  pj ,  that  is,  perform  the  first  order  approximation, 
we  have 


A  do  1  ^-4 

Apk  =  -pk - 2^ 

«0  Ofc-l  t^o 


k- 1 


Cfc— •  Cfc_,  2J  {^j—lJ—i^Pj  "b  C0jtj) 
;=i 


-  ePl*  (5.66) 
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If  we  do  not  count  the  ep,k,Za,j,i  and  Aao,we  end  up  with  the  (5.64). 

For  more  accuracy,  we  give  the  second  order  approximation  of  the  A  pk  as 


i  k— 1  k—1 

A  Pk  =  -  °k-i  H 

ak-l  ,= o  j=l 


3- 1 


A  pj  +  ea,j,i  +  Pj  ^  -i+«  A  pjl  +  ea,j1  ) 

ii=i 


Ao0  ea,i  1  ^-4  A 

—  Pk - Pk  - >  ajb-i,.-A  cfc_, 

«0  i=0  afc-1  ,=0 


(5.67) 


The  results  are  staggeringly  complex  for  continuing  this  procedure.  But  we  can  see  that 
the  errors  in  reflection  coefficients  at  kth  step  are  dependent  on  the  errors  in  previous  steps. 
Fortunatly,  we  do  not  need  to  pursue  the  higher  order  approximation  in  practice  since  the 
first  two  terms  dominate  the  errors. 


5.3.2  Finite  Analysis  of  Schur  Algorithm 

In  this  subsection,  we  are  going  to  continue  the  study  on  the  finite  analysis  for  the  Schur 
algorithm.  Schur  algorithm  is  given  in  section  5.1.2.  Since  the  erors  are  introduced  by  the 
finite  wordlength  of  computer,  the  algorithm  has  the  following  form  after  quantization 

hk-i,k 


Pk  = 


hk- 1,0 

hk,-j  =  hk-i  ,-j  +  pkhk-i,k+j  4-  eh,k,-j 

hk,k+j+l  =  Pkhk-\,-j-\  +  hk-i,k+j+l  +  Gh,k,k+j+l 

Taking  the  same  procedure  as  in  the  Levinson  algorithm,  we  have 

A  hk-i'k  +  P*A  hk-i,o 

A  Pk  - - 7 - eP,k 

tlk-1,0 

A  hk,~j  =  A  hk-i-j  +  hk~i,k+jA  Pk  +  Pk^k-i,k+j  +  ek-\,k+j 
A  hk,k+j+i  =  A  hk~i,k+j+ 1  +  ^t-i.-j-iA  pk  +  PfcA 


(5.68) 

(5.69) 

(5.70) 


(5.71) 

(5.72) 

(5.73) 


86 


We  expand  (5.72)  and  (5.73)  continuously  until 

k 

A  hk,~j  =  A  /tp.-j  +  5"!  (fe«-i.»+jA  p,  +  p«A  ~l~  eh,j,-j )  (5-74) 

«=i 

A  =  A0,fc+j+i 

k 

+  (^_i,_fc_ife+,_i  A  p,  +  pi  A  i,— j— *+,— i  +  Cfc,,-,fc+j+i )  (5.75) 

«=i 

is  reached. 

Let  us  compute  Ahk-i,o-  Setting  j  =  0  in  (5.69)  and  in  (5.13)  we  have 


Ajk.O  =  hk-1,0  +  Pkhk-i,k  +  eh,k,  o 

hk,  o  =  hk-i,o  +  Pkhk-i,k 


The  error  in  hk,o  is  derived  as 


A  hk, 0  =  A  hk-1,0  +  Pk A  hk-l,k  +  hk-l,k&  Pk  +  eh,k, o 


From  (5.71),  we  have  Ahk-i,k-  Substituting  it  into  Ahk,o,  we  get 

A  =  (l  -  p2)  A  /ijfe-i.o  -  Pkhk-i,o&  Pk  ~  Pkhk-i,o^P,k  +  hk-i,kA pk  +  eh,k, o 

Repeating  the  formula  again  on  A  hk- i,o,  A  hk-2,0,  •  •  •,  we  finally  have 

A  hk,o  =  A  h0,o  +  22  hk<o  f-2p,  A  p,  -  p,ePi,  +  (5.76) 

ho,o  i-o  V  / 


So  the  error  A  pk  for  Schur  algorithm  is 

APk  =  -p-A  Vo  +  pu  £  (2„,  A  A  +  ftep,)  -  pk  £  ^  -  e,.*  (5-77) 

ho,o 


t=0 


^».0  ^fc-1,0 
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Since  the  formula  is  very  complicated,  we  are  going  to  simplify  it.  Firstly,  we  note, 
from(5.74)  and  (5.75),  that  Ahk,-j  and  Ahk,k+j+i  have  the  form 

k 

A  hk,-j  =  ^2  fi{pji  PjPji  >  •  •  ‘)A  pj 
j=i 
k 

A  hk'lc+j+l  —  f zip  11  PlPji  >  •  •  ')A  Pi 

J=1 

If  we  only  consider  the  terms  that  are  not  dependent  on  Pj,  pjPj^ ,  •  •  we  have 


k-i 


A  hk-i,k  =  A  hotk  +  ^2  (hi~i,~k+iA  pi  +  &h,i,k) 

«=i 

Then  subsituting  Ahk-i^  into  (5.77)  ,  we  get  the  first  order  approximation  of  Apt  as 


A  pk  = 


hk- i,o  L 


k-l 

A  fto.fc  +  ^ 2  (hi-l,-k+iA  pi  +  £h,i,k) 

i=l 


eP,k 


(5.78) 


This  formula  is  much  simpler  than  the  one  in  Levinson  algorithm  (5.64). 

We  should  note  that  A  ho,*  =  Ac*  and  hfe_i,o  =  ctjfc-i-  How  about  ^“j1  A p,? 

Firstly,  from  (5.11),  we  write 


i— 1 

hi-l,-k+i  =  5 2  cj+k-iai-\,j  . 

i= o 


So 


fc-i  fc-i  i-i 

T~)  hj-i.—k+jA  pj  =  Cj+k-jOi-tjA  pj 

1=1  i=l  j=0 

k— 1  i 

=  5 2  Ck-jUi-lj-jA  pi 

«=i  i=i 


fc-l  Ac— 1 

=  ]£  C*-  H  A 

i=l  j=l 

This  result  has  shown  that  the  coefficient  of  Ap,  in  the  first  order  approximation  of 
A  pk  of  Levinson  and  Schur  algorithm  are  equal. 
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For  more  accuracy  we  can  include  the  terms  of  pi  as  the  second  order  approximation. 
From  (5.74)  and  (5.75),  we  have 

A  /ifc-i.fc  =  A  hd'k 


k- 1 


+  E 


t=i 


»-i 


hi-i'-j+k+i^iA  pi  +  pi  ho-k+i  +  52  ^H-i.ij-i+itApi,  +  eh,i,k-ij  +  en,i,k 


Substituting  A  A*-!,*  into  (5.77),  we  have 


1  k~1 
A  pk  =  t - 52  Pi 

"*-l,i  ,=0 


«-l 


A  C^_,'  +  )  — l,ij  — t'+fc  A  pix  +  Gh,ii,k—t 


Pk 

ho,o 


A  h0i0 


l  A  Pk  K  1  V'  C^|*.0 

+  A  /)fc,0  —  7 - A  /lo,0  —  pk\.  -7 - 

Ao.o  1=0  A.-.0 

where  A  pjti0  is  the  first  order  approximation  of  A  pk  given  in  (5.78). 
Now  we  will  prove  that 


k~\ 


^h,t,  0 


(5.79) 


Jl 


k-1  i-1  Ar-1  j- 1 

52 Pi  52  hi1~i<i1-i+kApil  =  52 pi  52(ah-'k*ii-i+i^Pi\) 

1=1  *i=i  *=o  ji=i 

Let  ji-j+i  =  u,  then 

*-i  k-i  i-i  k-i  i-i  fc-i+ji -j 

52 Ck-*52 Pi  52(ah-i&-i+i*ph) =  52 Pi  52  52  aii-i.«c^-«+ii-i^v. 

*=i  i=i  ii =1  i=i  ii=i  u=ii  -j 

Since  aj,_ 1)U  =  0  if  tt  >  ji  —  1  or  u  <  0,  we  have 


fc-i+ii  — i  ij  -i 

^  ^  ®ii—  i,»^fc— ■ »+ii  —  i  =  y  ]  aii-i,uCfc-u+ii  -i 
«=ii  -i  u=0 

Using  (5.11),  we  proved  the  equality. 

From  previous  study,  we  found  that  the  cofficients  of  A pi  of  the  errors  for  Schur  and 
Levinson  algorithms  are  same.  This  property  can  be  proved  by  induction.  This  also  means 
that  if  we  presume  c\  is  correct,  i.e.  A c±  —  0,  then  the  variance  of  errors  Apk  are  same 
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for  both  Schur  and  Levinson  algorithm.  Since  the  coefficiets  of  Ac,-  are  different  in  the 
two  algorithms,  in  general,  we  can  not  conclude  that  both  the  algorithms  provide  same 


accuracy. 


5.3.3  Finite  Analysis  of  Lattice  Algorithm 

For  the  lattice  algorithm,  we  take  the  same  procedures  as  in  previous  sections.  The  lattice 
algorithm  after  quantization  becomes 


n+fc—2 


9k  —  )  fk-l,ibk-U-l  "1*  Og,k 

,’= 1 

(5.80) 

Pit 

Pit  =  -  Opk 

Olk-\ 

(5.81) 

fk,i  =  fk-U  +  Pkbk-l,i-l  +  Oftk,i 

(5.82) 

bk,i  —  Pkfk-l,i  +  bk-l,i-l  +  Ob,k,i 

(5.83) 

at  =  ajt_i  +  9kpk  +  e0,k 

(5.84) 

The  errors  in  the  lattice  algorithm  are 


n+fc— 2 

y.  +  fk-l,iAbk-l,i-l )  +  Cg, k 

t=l 

AfiTfc  +  Pk£&k-l 

ep,k 


Agk 
Apk 
A fk,i  ■ 

&bk,i  ■ 

Aafc  =  Ac*fc_i  +  /jfcAgfjfe  +  p*  Ap*e0)* 


Olk- 1 
k 

A fo,i  +  y(^-i,i-iA°j  +  pjAAj-i,,--!  +  e/jt,- 

i=i 

k 

Abo,,-k  +  +  pjA/i-i.t-fc+i  +  e6j,«-fc+j 

i- 1 


(5.85) 

(5.86) 

(5.87) 

(5.88) 

(5.89) 
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Solving  for  Agk  in  (5.86),  substituting  it  in  (5.89)  and  simplifing  the  result  similarly  as  in 
previous  case,  we  have 

Aak  =  —A a0  4-  ak(-2piApi  +  —  +  p,eP).)  (5.90) 

<*0  -=o  °i 

At  this  point,  we  are  going  to  discuss  the  term  Agk  in  these  formula.  Firstly,  let  us  write 
the  Agi  as 

n— 1 

A^i  =  r(6o,.-iA/o,t  +  /o,:A5o,t-i)  +  eg>k 

1=1 

This  term  is  independent  on  Ap\. 

For  Ag?,  we  have 

n 

A^2  =  2(6i,i-iAfi,<  +  /i,«A5i,.-i)  +  e3i2 

«=i 

Using  (5.87)  and  (5.88),  we  have 

A/i,j  =  A/0,j  +  fto,i_iA/>i  +  piA5o,i-i  +  e/ti, ,• 

A6i,,-_i  =  Aio,«-2  +  /o,«-iApi  4-  piA/o,,_i  4-  e&, i,,-_i 

Using  the  lattice  algorithm  to  replace  the  terms  6i,j_i  and  /1,,  in  the  A^t,  we  have  the 
coefficient  of  Api  as 

n— 1 

[(Pl/o,«-l  +  &0,»-2)V«-l  +  (/o,«  +  Pi &0,t-l)/o,«-l] 

t=l 

Since  &o,«  =  fo,i  =  it,-,  <4  =  ^  J2]=o  UjUj-i  and  pi  =  —  we  obtained  a  zero  cofficient  of  Api. 

The  coefficients  of  Ap,  in  A^j.  are  zero.  This  is  am  advantage  of  lattice  algorithm  over 
both  Levinson  and  Schur  algorthms.  Keeping  this  in  mind,  we  are  not  going  to  consider 
the  term  Api  in  the  errors. 
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In  a  manner  similar  to  the  Levinson  and  Schur  algorithms,  the  error  can  also  be  classified 
by  terms  containing  products  of  pt.  The  first  order  approximations  are  given  as 

n+k— 2 


i  n-r*-, 

&Pk,  0  =  -  J2 

Qk- 1  j=l 


k 

I 

t'=l 


-f  ^2  +  /fc-l,j(AUj  +  E  eM,i-A:+«) 


i=l 


+  Cp.A:  (5.91) 


If  we  need  the  second  order  approximation,  that  is,  including  the  term  pi,  we  have 

„  fc-i  „ 

Aojt.1  =  ^Vfc.o - Ado  -  pk  2^  — 

ato  ~{  oa 


n+fc-2  f 


1  n+*— 

— — —  V 
"fc-i  L 


k—\ 


k- 1 


P,&bo,j-i  +  fk-\,j  p'Af°  J-k+i 


i=l 


i=l 


(5.92) 


From  this  result,  one  can  not  get  the  result  given  by  [12].  The  reason  is  that  we  removed 
the  term  A^*  in  (5.89)  by  (5.86).  If  we  removed  the  term  Apk  in  (5.89)  by  (5.86),  we  will 
get 


Aa*  =  (1  +  pD&cck- i  —  Ipk&Jk  ~  Pk<*k~iep,k  +  eQ,jt 


Ignoring  all  errors  except  of  Aafc_i,  we  have 

k—  1 

Aar*  =  Aax  JJ(1  +  p ]). 

j= 2 

Substituting  this  into  Ap*  of  (5.86),  we  have 

A Pk  =  — —Ac*!  H(i  +  Pj) 

Qk-i  j=0 

This  is  the  result  given  by  [12]. 

Comparing  the  results,  one  can  conclude  that  our  results  are  more  accurate  than  the 
one  given  in  [12].  Another  advantage  of  our  method  is  that  our  method  can  be  used  to  get 
more  accurate  results  by  including  the  terms  like  pip^ ,  and  so  on. 
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5.4  Finite  Analysis  of  Split  Algorithm 


5.4.1  Finite  Analysis  of  Split  Levinson  Algorithm 

The  errors  analysis  for  split  Levinson  algorithm  is  almost  same  as  those  for  Levinson  algo¬ 
rithm.  The  split  Levinson  algorithm  after  quantization  becomes 


“Tfe  —  ^  ^  Ci&k,i  4"  ^T,fc 

t=0 

0k  ~ 


(5.93) 

(5.94) 


—  Sk,i  +  —  0k^k-l,i-l 


(5.95) 


pk  =  1  — 


1  +  Pk- 1 


(5.96) 


The  errors  after  quantization  can  be  written  as 


A 7-fc  =  [Cj  As*,.-  +  s*,,-Acj]  +  eT,k 

i= 0 

**  Arfc  -  0k &Tk-i  , 

&0k  =  - +  e0,k 

Tk- 1 


(5.97) 

(5.98) 


As*,,'  —  +  Asfc_ii$_i  —  0k-l&Sk-7,i-\  +  •Sfc-2,«'-l  A^Jfc-1  +  es,k,i  (5.99) 


If  we  repeat  the  formula  As*,,-  to  Asfc_i,j,  As*.!,;.!,  •  •  by  induction,  we  can  get 

fc-l  k—j 

^k,i  =  (—0j&Sj-l,i-m  —  Sj-iti-m&0j  +  eSij+i,,_m+i)  (5.100) 

i=l  m=l 


5.4.2  Finite  Analysis  of  Split  Schur  Algorithm 


The  Schur  algorithm  after  quantization  has  the  form 


0k  = 

Ufe-1,0 


(5.101) 
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h  =  1 


A 


1  +  Pk-l 

Vk+ij  =  vk,j  +  6fc,j+ 1  -  $kVk-i,j+i 


The  errors  of  the  algorithm  after  quantization  are 


A?  : 

&Pk  =  (1  —  Pk) 

AiUfc.i  ■ 


Avk,o  -  /3kAvk-i,o 


vk-l,0 
A °k-i 


+  e/3  ,fc 


A3fc 


+  eP,fc 


1  +  Pk-l  1  +  Pk-l 
k— 1  fc— j  * 

E  E  err,  i,  (  Vj— l,t+m  A^j  +  ev,j+l,i+m— 1 ) 

j=l  m=l 

fe— 1 

+  (Aui,j+m-i  4-  Ayiii+m) 


(5.102) 

(5.103) 


(5.104) 

(5.105) 


(5.106) 


m=l 


5.4.3  Finite  Analysis  of  Split  Lattice  Algorithm 

The  lattice  algorithm  after  quantization  becomes 


2fk 

h 

Pk 


n+fc— 1 


«=0 


+  eT,fc 


1 -  +  eP,k 

Tk- 1 


1  - 


A 


1  4-  Pk-i 


+  eP,k 


tik,i  +  1  -  PkWk-i,i-i  +  eWtk+i,i 


The  errors  of  lattice  algorithm  are 


At* 


A % 

&Pk 


n+k—l 

wk,i&Wk,i  +  eT>k 

i=0 

Arfc  -  /3k&rk-i  , 

-  +  C<3,* 


Tfc-i 


(1  -  Pk) 


A>* 


A3fc 


1  +  /?*-!  1  +  /?fc+l 


+  eM 


(5.107) 

(5.108) 

(5.109) 

(5.110) 


(5.111) 

(5.112) 

(5.113) 
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(5.114) 


Aiojt.j 


fc- 1  k—j 

Z  E  CT-lii  -  wj.u-.A9f  +  e  tu,j+l,«— m+1 ) 

j=l  m=l 
fc— 1 

+  C'k'S1  (Au;,_m+i  +  Aivi-m) 

m=l 


Taking  the  same  steps  as  lattice  algorithm,  we  can  prove  that  At*,  is  independent  on 
Adj,  1  <  J  >  A:  —  1.  So  we  simplify  the  Atn*;,,-  into 


AtOfc,,- 


fc-i  fc-j 

E  E  c?-l  1  (  -f"  Cw,j+l,i— m+1 ) 

j= 1  m=l 


fc— 1 

+  C™- 2  (A^l.i-m+l  +  An;* ,t-m) 

m=l 


(5.115) 
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Chapter  6 


On  Reduced  Polynomial  Based  Split 
Algorithms 


One  of  the  fundamental  problems  in  linear  prediction  theory  is  to  compute  the  predictor 

polynomial  and  the  reflection  coefficients  (RCs)  for  a  given  real,  symmetric,  positive-definite 

Toeplitz  matrix  (henceforth,  referred  to  simply  as  Toeplitz  matrix).  For  such  a  task,  in  the 

context  of  order- recursive  algorithms,  the  Levinson-Durbin  algorithm  [1,  28]  was  the  most 

! 

computationally  efficient  algorithm  until  recently.  Recently,  an  algorithm  termed  ‘  the 
split  Levinson  algorithm  ’  has  been  described  in  [8],  which  requires  approximately  half 
the  number  of  MULT  and  same  number  of  ADD  as  the  Levinson-Durbin  algorithm.  The 
split  Levinson  algorithm  is  based  on  the  correspondence  between  the  set  of  polynomials 
orthogonal  on  the  unit  circle  and  the  set  of  polynomials  orthogonal  on  the  interval  [-1,  1) 
of  the  real  line.  It  is  also  closely  related  to  the  algorithms  for  zero  location  with  respect  to 
th  unit  circle  described  in  [27]. 
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The  Schur,  lattice,  and  the  normalized  lattice  algorithms  may  also  be  employed  for  the 
computation  of  the  RCs  for  a  given  Toeplitz  matrix  [29,  21,  30].  Each  of  these  algorithms 
possesses  its  own  set  of  desirable  features  and  is  closely  related  to  the  Levinson-Durbin 
algorithm.  In  view  of  such  a  close  relationship  and  the  existence  of  the  split  Levinson 
algorithm,  the  split  Schur,  lattice,  and  the  normalized  lattice  algorithms  have  also  been 
described  in  a  recent  paper  [9],  which  require  significantly  less  MULT  than  the  Schur, 
lattice  and  the  normalized  lattice  algorithms  respectively,  while  the  number  of  ADD  remain 
approximately  the  same. 

In  this  chapter,  we  derive  another  form  of  computationlly  efficient  three-term  recurrence 
relation  for  the  computation  of  RCs  (  and  the  predictor  polynomial  for  the  case  of  Levinson 
algorithm  )  associated  with  a  given  Toeplitz  matrix.  Following  the  terminology  introduced 
in  [8],  we  term  the  algorithms  as  the  split  algorithms.  These  new  split  algorithms  are  similar 
to  the  split  algorithms  of  [8,  9],  even  though  the  computations  involved  are  entirely  different. 
The  relationships  between  the  various  algorithms  are  also  established  in  the  paper. 

It  is  worthwhile  to  mention  here  that  the  three-term  recurrences  termed  the  split  Levin¬ 
son  algorithm,  derived  in  this  paper  and  in  [8]  were  perhaps  derived  in  [40]  for  the  first  time 
in  the  context  of  the  one-dimensional  inverse  problem  of  reflection  seismology.  However, 
their  application  to  computing  the  reflection  coefficients  and  the  Levinson  polynomial  in  a 
computationally  efficient  manner  was  demonstrated  in  [8]  and  is  being  done  here  as  well. 
Also,  the  present  work  and  [8]  extend  these  computational  savings  to  the  Schur,  lattice, 
and  normalized  lattice  algorithms.  Finally,  in  this  paper,  we  also  describe  a  new  lattice 
realization  of  a  digital  filter  and  a  new  computationally  efficient  algorithm  for  solving  the 
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Toeplitz  linear  system  Tx  =  b. 

This  Chapter  is  organized  as  follows.  Section  6.1  sets  up  the  notation  for  the  sequel 
and  develops  the  necessary  mathematical  background  for  the  work.  In  Section  6.2,  we 
introduce  the  new  split  Levinson  algorithm  and  present  computationally  efficient  algorithms 
for  the  computation  of  the  RCs  and  the  predictor  polynomial.  The  various  interrelationships 
between  the  new  split  Levinson  algorithm  and  the  previous  Levinson- Durbin  and  split 
algorithms  are  also  established.  In  Section  6.3,  the  corresponding  split  Schur  algorithm 
is  derived,  and  in  Section  6.4,  we  derive  the  corresponding  lattice  and  the  normalized 
lattice  algorithms.  These  algorithms  lead  to  a  new  lattice  structure  for  finite-impulse-re- 
sponse(FIR)  filtering  of  discrete  data.  Based  on  the  new  split  Levinson  algorithm,  a  feist 
algorithm  for  solving  the  general  linear  system  Tx  —  b  is  briefly  outlined  in  Section  6.5. 
It  is  worthwhile  to  mention  here  that  in  order  to  avoid  introducing  new  mathematical 
terminology,  a  simple  matrix-vector  approach  has  been  adopted  to  describe  the  various 
algorithms. 


6.1  Notation  and  Classical  Algorithms 

The  mathematical  preliminaries  fundamental  to  the  content  of  this  work  are  presented  in 
this  section.  The  various  notations  are  defined,  and  certain  well  known  results  on  Levin¬ 
son,  Schur,  lattice,  and  the  normalized  lattice  algorithms  and  their  connection  to  Toeplitz 
matrices  are  described  [1,  28,  29,  21,  30]. 
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6.1.1  Notation 


Let  1Z  denote  the  field  of  real  numbers  and  T  be  a  Toeplitz  matrix  of  order  (n+1),  specified 
in  term  of  its  elements  in  the  first  column,  i.e.,  c,e7?.,  i  =  0,1,  •••,«.  For  a  given  matrix 
T,  let  Tk  denote  the  principal  submatrix  of  T  of  order  (k+1).  Note  that  Tn  =  T.  A 
(k+l)-dimensional  cloumn  vector  x*  =  (xk,oXk,i  ••’xk,k)i  can  equivalently  be  expressed  as 
a  polynomial  Xk(z)  =  Yli=oxk,iZ'  of  degree  k  and  vice  versa.  The  equivalence  between  the 
vector  notation  x*  and  polynomial  notation  Xk{z)  is  used  throughout  this  work  in  order  to 
simplify  the  form  of  the  expressions.  In  context  of  filtering,  z_1  may  be  considered  as  the 
unit  delay.  Let  J  denote  the  exchange  matrix  having  ones  on  the  antidiagonal  and  zeros 
elsewhere.  The  order  of  the  matrix  J  is  understood  to  be  conformable  with  the  context 
in  which  it  appears.  For  a  given  polynomial,  ijt(z)  =  xk,iz\  let  x*(z)  denote  the 
reciprocal  polynomial,  x*(z)  —  Yli=oxk,k-iZl-  In  vector  notation,  we  have  x*  =  Jxk.  A 
vector  Xfc  (polynomial  xjt(z))  is  called  symmetric  if  x*  =  Xjt  (xjt(z)  =  X/t(z))  and  anti¬ 
symmetric  if  Xfc  =  —x*  (xfc(z)  =  —  xjt(z)).  In  either  of  the  two  cases,  x*  is  completely 
determined  by  approximately  half  the  elements  (|  if  k  is  even  and  (^r-)  if  k  is  odd).  For 
a  given  vector  x*  =  ( Xk<0  •  ■  •  x*,*)*, vector  xm  =  (x*,,-  •  •  •  x*,, •+„,)'  is  called  a  subvector  of  x*. 
Of  course,  i  >  0,  and  i  +  m  <  k.  For  i  =  0,  xm  is  the  principal  subvector  of  x*,  and 
for  i  =  (k  —  m)/ 2,  xm  is  center-subvector  of  Xk  (k-m  is  even  ).  Any  vector  x*  can  be 
written  as  sum  of  a  SYM  vector  s k  and  an  ASYM  vector  a*,  where  s*  =  \[xk  +  X*]  and 
a*  =  i(xfe  —  Xfc]  =  x*  —  s*.  If  xm  is  the  center-subvector  of  x*,  then  xm  is  the  center- 
subvector  of  Xfc.  From  this  discussion,  it  is  clear  that  center-subvector  of  x*,  i.e.,  xm, 
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can  be  written  as  sum  of  a  SYM  vector  sm  and  an  ASYM  vector  am,  where  sm  and  am 
are  the  center-subvectors  of  s&  and  a^  respectively.  We  use  Ylx  t°  denote  the  sum  of  the 
elements  of  x*,  i.e.,  Ylx  —  Yli=oxk,i,  and  ||x;t||  is  used  to  denote  Euclidean  norm  of  Xjt, 
i.e.,  ||xfc||  =  (J2i=oxl,i)2  ■  The  Toeplitz  matrix  Tk  satisfies  the  properties  1)  Tk  =  T£,  2) 
JTkJ  =  Tk. 


6.1.2  The  Levinson-Durbin  Algorithm 

The  predictor  polynomials  <f>k(z)  satisfy  the  two-term  recurrence  relation, 

<j>k(z)  =  <j>k-i{z)  +  Pk*ik-i{z)  (6.1) 


where 


and 


fc— l 

Pk  =  —Vk- 1  H  ck-i<t>k-\,i 

iss  0 


(6.2) 


cr k  =  cr k-i(l  -  pi)  k  =  l,2,---,n 


(6.3) 


Here,  at  =  ,  cr~\  =  1,  and  <j>o(z)  =  1.  The  scalars  pi,/>2,  •  •  • ,  pn  are  called  the  RCs.  It 

is  clear  from  (  6.1)  that  (f>k,o  =  4>k-i,o  =  •••=*  <t> o,o  =  1,  and  <f>k,k  =  Pk-  The  corresponding 
vector  4>k  satisfies  the  system: 


Tk(f>k  =  [crk  0  •  •  •  0  ]* 


(6.4) 


and,  from  (  6.2),  d’fc-i  satisfies  the  system: 


Tk 


<t>k- i 

0 


=  [<7fc_i  0  ■  •  •  0  -  pkCTk-lf 


(6.5) 
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Also,  substituting  z=l  in  (  6.1),  we  get, 


<M1)  =  ^fc-i(l)  +  Pk<i>k-i{l) 

Since  <£fc-i(l)  =  we  have,  <£fc(l)  =  (1  +  pk)<t>k-i{l),  or 

wi)=n(i+«>  <6-6) 

*=i 

The  above  expression  was  first  obtained  in  [6].  A  straightforward  analysis  shows  that  this 
algorithm  requires  n2  -+•  0(n)  MULT  and  n2  4-  0(n)  ADD. 


6.1.3  The  Schur  Algorithm 

In  essence,  the  Schur  algorithm  performs  the  triangular  decomposition  of  the  Toeplitz 
matrix.  The  Schur  polynomials  (SPs)  satisfy  the  two-term  recurrence  relations, 


z0n_fe_i(z)  =  9n-k(z)  +  Pk^n-k(z) 


(6.7) 


A„_fc_i(z)  =  PkOn-k(z)  +  A n-k(z) 


(6.8) 


where 


9n-k,  0  t.  i  n 

Pk —  a  k  —  1^2, 


(6.9) 


An_/t,o 

The  initial  SPs  0„_i (z)  and  An_i(z)  are  specified  as  9n-\{z )  =  c<+ \z\  and  A„_i(z)  = 

££?  c>z'-  The  connection  between  the  Levinson-Durbin  and  the  Schur  algorithm  is  as 
follows.  For  the  linear  transformations, 


u?>  =  rw_i  0  •  •  •  0)' 


(6.10) 
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and 


v?>  =  r(0 0 ...  0)'  (6.11) 

the  SPs  «„-k(z)  and  A„-t(z)  are  defined  as  9„_t(z)  =  £"T0‘  »S+,z',  and  A„_t(z)  = 
4_/,=o  vn,fc+»2  •  Note  that  the  superscript  (k)  has  been  used  here  to  demonstrate  the  de¬ 
pendence  of  the  various  quantities  on  the  index  k.  Using  (6.1),  (6.4)  and  (6.5),  we  obtain 
the  recurrence  relations  in  (  6.7)-  (6.9).  Using  (6.4)  and  (6.5),  it  is  easy  to  see  that  we  also 
have  the  identities, 

u!,0  =  uiJ(An-fc,o)  =  (Tk-I 

«i*iw «-».o)  =  «S3  =  —pk&k~i 

and 

«S  =  ^i=0,  *  =  1,2,  •  •  •  —  1 

The  complexity  of  the  Schur  algorithm  is  same  as  that  of  the  Levinson-Durbin  algorithm. 


6.1.4  The  Lattice  Algorithm 


In  certain  signal  processing  applications,  the  Toeplitz  matrix  appears  as  the  product  T  = 

AA\  where  A  is  an  (n  +  1)  x  (N  +  n)  order  upper  triangular  Toeplitz  matrix  of  rank  n+1, 
given  by, 


A  = 


®o  «i  an  •  •  •  o^/_!  0 

0  a0  •  •  •  an_j  •  •  •  ajv-2  ajv-i 


0 

0 


(6.12) 


0  0 


Qo  *  •  ■  a/v_n_  i  as-n  •  •  •  a;v- 1 
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Here  N  >  n.  Let  Ak  be  the  (k  -f  1)  x  (jV  -f  k)  order  principal  submatrix  of  A,  then  we  have 
7*  =  AkA‘k,  k  =  0,  •  •  •  ,n.  The  lattice  polynomials  (LAPs)  satisfy  the  two-term  recurrence 
relation, 

CN+k-l(z)  =  ev+k-liz)  +  ZpkfN+k-2(z)  (6.13) 

and 

fN+k-i(z)  =  pkeiv+k-2(z)  +  zfN+k-2{z)  (6.14) 


where 

Pk  =  _ °rJfe-l(^+fc-20) 
In  addition,  we  also  have 


'  0  X 


^  e^+k-2  ) 


k  =  1, 2,  •  •  • ,  n 


(6.15) 


°k~  i  =  ||fiv+fc-2|l  =  ||ejv+fc-2||  (6.16) 

The  initial  LAPs  are  specified  as  ew_i (z)  =  / n-i{z )  =  a»^,•  The  connection  between 

the  Levinson-Durbin  and  the  lattice  algorithm  is  as  follows.  For  the  linear  transformation, 


Ujv+t-i  =  Kktfk- 1 


(6.17) 


and 

v^+*- i  =  4H(0 


the  LAPs  e^+t_2(z)  and  fs+k-2{z)  are  defined  as 

N+k- 2 

e^+k-iiz)  =  ^2  UH+k-uz\ 

i=0 


(6.18) 


and 


N+k— 2 

fN+k-2(z)  =  ^  Vjv+fc-i,  ,+iz'. 

t=0 


103 


Using  (6.1),  (6.4)  and  (6.5),  we  obtain  the  relations  in  (6. 13)— (6. 16) .  Also,  premultiplying 
both  sides  of  (6.17)  by  Ak  ,  we  get, 

AfcUjv+fc-i  =  Tk(<f>k-i  0)‘  =  0  •••  0  (6.19) 

As  a  result,  the  computation  of  pk  in  (6.15)  may  alternatively  be  performed  as, 

N-2 

Pk  =  —<Tk-\  aieH+k-2,i+k  (6.20) 

is  0 

However,  the  expression  in  (6.15)  is  used  more  frequently.  Computation  of  RCs  using  the 
lattice  algorithm  in  (6.13)-(6.15),  requires  3 JVn  +  1.5n2  +  0(N  +  n)  MULT  and  ADD.  The 
recurrence  relations  may  also  be  used  to  obtain  the  lattice  realization  of  the  FIR  filter 
having  transfer  function  (f>(z~l).  Such  a  realization  is  shown  in  Fig.l  for  order  n=5. 


6.1.5  The  Normalized  Lattice  Algorithm 


The  normalized  lattice  algorithm  is  the  normalized  version  of  (6.13)-(6.15),  where  all  the 
vectors  have  unit  norm.  Using  (6.16),  the  normalized  lattice  algorithm  may  be  written  as 


9N+k- i(z)  =  (1  —  Pk)  ^{9N+k- 2(2)  +  ZPkhN+k-liz)} 


hf/+k-l(z)  =  (1  -  pl)~*  [Pk9N+k-2(z)  +  zhtf+k- 2(2)] 


(6.21) 

(6.22) 


where 


Pk  =  (gN+k-2  0) 


0 

hjV+Jfc-2 


(6.23) 


Here,  g^+*_2  =  jjggggfl  and  h^+*- 2  =  ||fR»Id| >  and  the  initialization  is  given  by  gN-i(z )  = 


hrt-i(z)  —  (Silo1  ai)~*  Silo1  aiz'-  The  normalized  lattice  algorithm  requires  5iVn-f2.5n2+ 
0(N  +  n)  MULT,  3 Nn  +  1.5n2  +  0(N  +  n )  ADD,  find  n  square-root  computations. 
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We  conclude  this  section  by  stating  that  the  computationally  efficient  alternatives  to 
the  above  algorithms,  which  are  presented  in  the  following  sections,  follow  a  very  similar 
approach  for  the  various  computation.  The  main  difference,  perhaps,  is  in  the  use  of  a 
three-term  recurrence  in  place  of  a  two-term  recurrence. 


6.2  A  New  Split  Levinson  Algorithm 


In  this  section,  we  present  a  new  split  Levinson  algorithm  based  on  the  recursive  evaluation 
of  a  symmetric  polynomial  using  a  three-term  recurrence  relation.  Consider  the  system  of 
linear  equations: 

TkC I*  =  [Ok  *  •  •  A*]4  (6.24) 


Where  the  scalar  Ok  is  chosen  such  that  qk,o  —  1*  Clearly,  q*  is  symmetric ,  and  ,  therefore, 
qk,k  =  1.  The  above  system,  though  relatively  unknown  in  signal  processing  applications, 
appears  frequently  in  problems  of  discrete  inverse  scattering  [31,  32].  In  the  present  context, 
we  are  interested  only  in  an  efficient  computation  of  the  RCs  and  the  predictor  polynomial. 
To  obtain  a  recurrence  for  qk(z )  we  proceed  as  follows.  Writing  (6.24)  in  terms  of  k-1  and 
k-2,  it  is  straightforward  to  see  that, 


Tk 


Tk 


qfc-i 

o 

o 


qjfe-i 


=  [#*-i  •  •  *  0fc- 1  Tfc-i]* 


=  [7fc-i  9k- 1  •  •  •  fljfc-i]* 


(6.25) 


(6.26) 


105 


and 


Tk 


0 

q*-2 

0 


=  hk-2  Ok-2  •  •  •  Ok-2  7fc-a]* 


Based  on  (6.25)-(6.27),  the  recurrence  for  q*  can  be  written  as, 


q* 


“ 

■ 

q*:-l 

+ 

0 

0 

qjb-i 

o 

qfc— 2 

o 


Premultiplying  both  sides  of  (6.28)  by  Tk,  we  get, 


0k 

0k-l 

7Jt-i 

7fc-2 

0k 

0k- 1 

+ 

0k-\ 

rH 

i 

3 

1 

0k- 2 

0k 

0k- 1 

0  k- 1 

0k- 2 

0k 

7*-i 

0k-\ 

lk-2 

The  above  equation  is  satisfied  if  the  scalar  ak  is  such  that 


0k- 1  +  7*-i  —  Qk7k-2  —  20fc_i  -  ak6k- 2 


(6.27) 


(6.28) 


(6.29) 


(6.30) 


or 


Ok- i  -  7fc-i 
0k-  2  ~  Ik-2 


Once,  ak  is  known,  0k  is  obtained  as, 


0k  =  2  0k-i  —  ak0k-2 


(6.31) 


(6.32) 
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Also,  (6.28)  can  be  written  in  the  polynomial  form  as, 


qk(z)  =  (1  +  z)qk-i(z)  ~  akzqk.2(z) 


(6.33) 


The  recursive  computation  of  qk(z )  from  qk-i(z),  k  =  2,  •••,«,  consists  in  (i)  computing 
7it_i  using  (6.25),  (ii)  computing  ak  using  (6.31),  (iii)  computing  qk{z)  using  (6.33),  and  (iv) 
computing  0k  using  (6.32).  The  initial  values  are  given  by  q0(z )  =  1,  0Q  =  Co,  qi(z )  =  1  +  2, 

=  co  +  ci,  and  71  =  c2  +  Ci.  Using  the  symmetry  property  of  qk(z),  qn(z )  can  be  « 
computed  recursively  in  0.5n2  +  0(n )  MULT  and  n2  +  0(n)  ADD.  In  the  following, we 
establish  the  relationships  between  the  polynomials  qk(z),  predictor  polynomials  <t>k(z),  and 
the  polynomials  used  in  the  split  Levinson  algorithm  of  [8].  We  consider  the  relationship 
between  qk(z)  and  <f>k(z)  first.  Combining  (6.24)  and  a  scaled  version  of  (6.26)  (scalar  = 


-Afc),  we  get, 


Setting 


Tk{  q*  —  A* 


0 

q*-i 


0*  —  Ajfe7fc_i 


0k  —  Afc0jt_i 


Ok  ~  Afc0fc_i 


(6.34) 


(6.35) 
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(6.34)  reduces  to, 


Tk{ 


qjt  —  Afc 


0 

q*-i 


>  = 


— rk-\)6k 
0*-l 

0 


(6.36) 


A  comparison  of  (6.4)  and  (6.36)  reveals  that, 


Mz)  =  <lk(z)  -  Xkzqk-i{z) 


(6.37) 


and 


0k 


(flfe-i  ~  lk-\)0k 
Ok- 1 


(6.38) 


Since  —  1  and  the  first  element  of  q*  in  (6.36)  is  also  1.  Also,  (j>k,k  =  pk  which,  in  turn, 
implies  that  (  refer  to  the  r.h.s.  of  (6.36)  again), 


</>k,k  =  Pk  =  1  —  A* 


(6.39) 


or. 


Xk  =  l  ~  Pk 


(6.40) 


Substituting  (6.35)  in  (6.38),  we  get, 


a  °k 

Ok- i  -  Ik-i  =  Yk 


Combining  (6.35)  and  (6.40),  it  can  be  stated  that, 


A 

0k- 1 


=  1  —  Pk 


(6.41) 


(6.42) 
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or,  in  a  non-recursive  form,  we  have, 


k 

Ok  =  Co  JJ(1  —  pi) 

1=1 

Dividing  both  sides  of  (6.32)  by  Ok-  1  and  simplifying  using  (6.42),  we  also  get 

ajfc  =  (1  —  Pfc-i)(l  +  Pk) 


(6.43) 


(6.44) 


Thus,  an  order-recursive  algorithm  for  the  computation  of  the  predictor  polynomial  4>n(z) 
and  the  RCs  pi,  />2,  •  •  • ,  pn  consists  in  computing  qk(z),  k  =  2,  •  •  • ,  n  recursively  as  outlined 
above,  computing  pk,  k  =  1,2,  •••,n  using  (6.39),  and  then  recovering  <t>n(z)  using  qn(z ) 
and  qn-\(z)  as  in  (6.37).  This  new  algorithm  is  summarized  in  Table  I.  The  recovery  of 
(pnyz)  requires  an  additional  0.5n  MULT  and  n  ADD.  Of  course,  if  one  is  interested  in 
computing  the  RCs  only,  then  (6.37)  may  not  be  used.  Also,  if  it  is  required  to  compute 
all  the  predictor  polynomials  ^fc(z),  k  =  1,2, •  ••  ,n,  then  the  procedure  descibed  here  may 
still  be  used.  However,  such  a  procedure  would  then  require  0.75n2  +  0(n )  MULT  and 
1.5n2  +  0(n)  ADD. 


Another  important  relation  is  obtained  by  setting  z=l  in  (6.37)  and  considering  the 
polynomials  q'k(z )  =  0^lqk{z).  This  leads  to 

^fc(l)  =  Ok  ^2  —^kOk-i 

q'  q' 

Using  (6.6),  (6.35),  and  (6.43),  we  get, 

k  k—1  k  i  I 

1  +  Pi 


E-E-fnfr* 

q'  q'  «=1  1  P* 


or 


q' 


(6.45) 
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as  q'0  =  Cq1.  The  above  relation  has  been  employed  in  problems  of  discrete  inverse  scattering 
[31,  32,  33]  and  in  problem  of  computing  likelihood  ratios  for  known  signals  in  stationary 
correlated  Gaussian  noise  [34].  Note  that  the  vector  is  solution  to  the  system, 

Tkqk  «  [1  1  1]*  (6.46) 


and 


I>(ii  •••  ill'll i  i]‘ 

q ‘ 


(6.47) 


Since  q'k(z)  =  6klqk(z),  we  also  have, 


k 

£  = 
q 


iic1  _ p*) + nc  (n(i + Pt)  n 

,*=i  j= i  i«=i  i=j+i  j. 


(6.48) 


Furthermore,  substituting  z=l  in  (6.33),  we  get  the  recursion, 

k  k-l  k- 2 

E  =  2  E  -ak  E 

q  q  q 


(6.49) 


We  now  turn  to  the  relationship  between  qk(z)  and  the  polynomials  used  in  the  split  Levin¬ 
son  algorithm  of  [8].  Adding  (6.25)  and  (6.26)  and  subtracting  a  scaled  version  of  (6.27) 
(taking  into  account  Ak_i  =  we  must  set  the  scaler  =  2Ajt— i )  we  get, 


r 

- 

qt-i 

0 

Tk* 

0 

+ 

q*-i 

—  2Afc_j 


0 

qjfe-2 

0 


}  =  frt  o  •  •  •  o  Tky 


(6.50) 


where 


Tk  =  Ok-i  +  7*-i  “ 


2^fc-i7fc-2 
Ok- 2 


(6.51) 


no 


Let, 


P*  = 


qfc-i 

c 


+ 


o 

q*-i 


—  2Ajt_x 


0 

qfc-2 

0 


In  polynomial  form, 


pk(z)  =  (1  +  z)qk_i(z)  -  2A*_i zqk-.2{z) 


Using  (6.30),  (6.38),  (6.42),  and  (6.44),  it  may  be  established  that, 


n  = 


Ok 


1  +  Pk 


and  (6.50)  can  be  written  as, 


(6.52) 


(6.53) 


(6.54) 


TkPk  =  [rife  0  •  •  •  0  t*]* 


(6.55) 


It  is  clear  that  the  vector  p*  is  symmetric,  and  pk,o  =  1.  A  comparison  of  (6.54)  and  (6.55) 
derived  above  with  (6.28)  of  [8]  reveals  that  the  polynomials  pk{z)  as  expressed  in  (6.53)  is 
same  as  the  symmetric  polynomial  pk(z)  of  [8]. 

Also,  subtracting  (6.26)  from  (6.25),  we  get, 


=  (*fc-i-7*-i)[10 


•••  0  -1]‘ 


(6.56) 


Let 


or,  in  polynomial  form, 


“  ■ 

• 

q*-i 

0 

0 

q*-i 

Pk(z)  =  (1  -z)qk-i{z) 


(6.57) 


(6.58) 
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and 


Tfc*  =  Bk- 1  -  7*- 1  (6.59) 

Combining  (6.41)  and  (6.59),  we  also  have 


Thus,  (6.56)  may  be  written  as, 


lip;  =  K  0  •  •  •  0  -r;]‘  (6.60) 

It  is  clear  from  (6.57)  that  the  vector  p *k  is  antisymmetric,  and  pkQ  =  —  p*k  k  =  1.  A 
comparison  of  (6.58)  and  (6.60)  derived  above  with  (6.21)— (6.23),  (6.30)  and  (6.31)  of  [8] 
reveals  that  the  polynomials  qk(z)  studied  here  are  the  same  as  the  reduced  polynomials 
qk(z )  of  [8].  The  analysis  performed  here  also  invalidates  some  of  the  statements  made  in 
[8]  regarding  the  reduced  polynomials. 

The  numerical  stability  of  the  algorithm  which  processes  qk(z)  recursively  in  place  of 
<j>k{z )  as  in  the  Levinson-Durbin  algorithm  is  an  important  issue  and  needs  further  analysis. 
Such  an  issue  goes  beyond  the  scope  of  this  paper.  Also,  it  is  worthwhile  to  mention 
here  that  the  recovery  of  predictor  polynomial  <f>n(z)  from  the  polynomials  studied  in  [8] 
requires  division  by  (1-z),  which  is  not  the  case  for  recovery  of  <t>n{z)  from  qn{z)  and  qn-\{z). 
Therefore,  it  may  be  preferable  to  use  the  algorithm  in  terms  of  qk{z),  k  =  1,2,  •  •  •  ,n. 

Using  a  similar  procedure  as  described  above,  we  can  obtain  a  three-term  recurrence  for 
the  polynomial  r*(z)  satisfying  the  system, 

TkTk  =  [0k  —Ok  •••  (-l)*0fc]‘  (6.61) 
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Note  that,  r*  =  (  — l)fcrjt,  i.e.,  f/t  is  symmetric  if  k  is  even  and  it  is  antisymmetric  if  k  is  odd. 
It  is  interesting  to  note  that  the  polynomials  rk(z)  may  be  obtained  from  the  polynomials 
qk(z)  via  the  change  of  variables  z  — *  —z.  The  polynomials  rk(z )  satisfy  the  following 
three-term  recurrence, 

rk{z)  =  (1  -  z)rk-i(z)  +  ackzrk-2(z) 

Finally,  we  conclude  this  section  by  stating  that  the  polynomial  qk{z)  is  a  monic  polynomial 
and  its  recursion  may  be  termed  the  ‘monic  recursion’.  By  an  appropriate  scaling  of  qk{z), 
one  may  also  obtain  the  ‘  balanced  recursion  ’  and  the  ‘  dual  recursion  ’[35].  Without  any 
loss  of  generality,  we  focus  on  the  monic  recursion  only. 


6.3  The  Split  Schur  Algorithm 

In  this  section,  based  on  the  relationships  between  the  classical  Levinson-Durbin  and  Schur 
algorithms,  we  derive  the  split  Schur  algorithm  corresponding  to  the  polynomials  qk{z)  for 
the  computation  of  the  RCs. 

Consider  the  linear  transformation, 

uW  =  T(q‘_10  •••  0)‘  (6.62) 

and 

vW  =  T( OqUO  •••  0)‘  (6.63) 

From  (6.25)  and  (6.26),  it  is  clear  that 

=  »&i  =  i  =  0,l,  •,*- 1  (6.64) 
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(6.65) 


„(*)  _  (*)  _ 

Un,k  ~  Vn, 0  -  7*-l 


uw  -  v(k) 

un,  i  —  un,»+l 


i  =  k,  •  •  •  ,n  —  1 


(6.66) 


Now,  we  define  the  split  Schur  polynomial  as  an_jt(z)  =  5Z£Zo'  un,l+izt-  As  a  result  an_jt,o  = 
fit-i-  Using  the  linear  transformation  and  the  various  identities  in  (6.62)-(6.66),  the  split 

Levinson  recurrence  in  (6.28)  is  transformed  into  the  split  Schur  algorithm.  Since  the 

* 

derivation  is  straightforward,  we  simply  summarize  the  various  expressions  in  the  following: 


=  Gin— it, o 

flfc-i  ~  7fc-i 
0Jt— 2  “  7fc-2 


zan-k-i(z)  =  (1  +  z)an.k(z)  ~  akan.k+i(z) 


(6.67) 


9k  —  20k-i  —  ctkOk-2 


Pk  =  1  -  2 -  k  =  2, 


The  initialization  is  given  by  a„_x(z)  =  £”=<)  ci+iz\  an_2(z)  =  £"_T0  (c,+i  +  c«+ 1)2'-  The 
above  described  split  Schur  algorithm  requires  0.5n2  +  O(n)  MULT  and  n2  +  0(n )  ADD.  It 
is  interesting  to  note  that  the  relationship  between  the  split  version  of  the  Levinson-Durbin 
and  Schur  algorithm  is  same  as  that  of  their  classical  versions. 


6.4  The  Split  Lattice  Algorithms 


In  this  section,  we  derive  the  split  lattice  and  the  split  normalized  lattice  algorithms  cor¬ 
responding  to  the  polynomials  g*(z)  for  the  computation  of  RCs.  Given  T  =  and 
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Tk  =  AkA‘k,  consider  the  linear  transformation, 


UAr+fc_i  =  i4*(q*-i  0)* 


(6.68) 


The  split  lattice  polynomials  (  SLAPs  )  are  defined  as,  fw+k-2(z)  =  2  uN+k-i,iz'  or, 

?N+k-2 


0 


=  u;v+fc_i.  From  the  above  transformation,  it  is  clear  that, 


0 

ffV+fc-2 


=  Ai 


0 

qA:-l 


Also,  since  Tk  =  AkAk,  premultiplying  both  sides  of  (6.68)  by  Ajt,  we  get, 


AfcUjv+fc-i  =  7fc(qJ._1  0)*  =  [0fc_i  •  •  •  Ok— i  TAr— i]4 


(6.69) 


(6.70) 


or 


N- 2 

7 k-1  =  5Z  aifN+k-2,i+k 

i=0 


(6.71) 


One  form  of  the  new  split  lattice  algorithm  is  obtained  by  using  the  linear  transforma¬ 
tions  in  (6.68)  and  (6.69)  on  (6.28).  Such  a  derivation  is  based  on  (6.71),  and  the  expressions 
for  the  corresponding  split  lattice  algorithm  can  be  summarized  as, 


N-2 

7*-l  =  51  QifN+k-2,i+k 

i=0 

0k- 1  -  7fc-i 

Qk  —  7 - 

0k- 2  ~  Ik-2 

fN+k- 1(«)  =  (1  +  z)ftf+k-2(z)  ~  <XkzfN+k- 3(z) 


0k  =  20k-i  —  ctkOk- 2 

0k 

Pk  =  1  ~  2 —  k  =  2,  ■  •  • ,  n 

0  k- 1 


(6.72) 
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The  initialization  is  given  by  /;v-i(-z)  =  E^o*  aiz\  In(z)  =  E£Lo(at+a.-i)z*,  9o  =  Eilo1  ab 
Another  form  of  the  split  lattice  algorithm  may  be  obtained  by  considering  the  inner  prod¬ 
uct, 


~  [fw+fc-2  0] 


qLi  ol  t„ 


0 


ftf+Jfe-2 


=  [q*-i  o]  AkA[ 


0 


qi-i 


qfc-i 


=  [q*-i  0  hk-i  Ok-i  ■  •  •  0/t-iF 


Since  (jk—i,o  =  1,  the  above  expression  simplifies  to, 


fc-i 


T/t  =  7fc_i  +  ~1) 

q 


or 


fc-i 


9k-i  —  7fc-i  =  0fc-i  5Z  ~^k 
q 

As  a  result,  the  split  lattice  algorithm  based  on  the  above  expression  may  be  written  as, 


Tfc  —  [f^+*_2  0] 


0 

^N+k-2 


9k- 1  —  7fc-i  = 
ak  = 
fN+k- l(z)  = 
9k 

k 

E 

q 


fc-i 

9k-i  ^  —Ik 
q 

9k- i  -  7fc-i 
9k-2  —  lk-2 

=  (1  +  z)fN+k-2(z)  -  akzfN+k.3(z ) 

=  2&k-i  -  Ctk0k- 2 

fc-l  k—2 

=  2j2~akJ2 

q  q 

9k  ,  0 

k  =  2,  ■  •  •  ,n 


(6.73) 


Pk  =  1- 


0*-x 
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The  initializations  are  £q  =  1,  =  2. 

The  normalized  split  lattice  algorithm  may  be  based  on  either  (72)  or  (73).  In  the 
following,  we  focus  on  (6.73).  Consider  the  norm, 


M  =  || — 2 1| 2  =  fN+k-2^N+k-2  —  [fw+fc-2  0] 


ftf+k-2 


Let  us  define, 


—  U;v+j._1U/V+fc_i  —  [Qfc_i  0]7fc 


=  [q£_x  O][0*-i  •••  &k- 1  7*-i]‘  =  0k- i 


f/V+fc-2  = 


?N+k- 2 
||  f/V+fc — 2 


Now  consider  the  inner  product, 


Tjfc-1  =  [fN+k-2  0] 


?N+k-2 


»>-i  si-1 -(»*-! 

h- 1  Ej-1 


Similarly, 


0k- 1  -  7*-i  ,  t 


0k- 2  ~  lk-2  ,  ~ 


Based  on  the  above  two  relations,  we  define  two  scalars  a^-x  and  bk-i  as 


_  I  -  Jfc-1  _  (0k- 1  -  lk-l)0k-2llq  2 
k  1  1— Jfc-2  (0k-  2  —  7ifc— 2)0Jfc— 1  Hq-1 


U-jEJ'v 
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Then  a*  =  ak-ib2k_l.  Replacing  each  of  f,-  by  |jfj|[f,  in  (6.73),  we  obtain  the  normalized  split 
lattice  algorithm.  The  various  expressions  for  such  an  algorithm  can  be  summarized  as, 


0 

lk-1  =  [f/V+fc— 2  0] 

■  _ 

fv+jfc- 2 

«*-i  =  (1  —  XJfc_1)(l  -  Jit-z)-1 

Oik  =  Ojt-l&fc-i 

0k  =  26k-l  —  Oik-20  k— 2 

k  it— 1  k-2 

£  =  2  £-<**£ 
q  q  q 

fc  *-i  1  ^ 

^  E(^-i  5Z)-1 

q  q 


(6.74) 


h 

fN+k-i(z)  =  +  z)fN+k- 2(z)  -  ajt-i 6jk_i /iv+Jt-3 ( z )] 

Pk 


0k 


1  _  =  2, •  •  •  ,n 

"t-i 


The  initializations  are  same  as  in  the  split  lattice  algorithm  except  that  the  initial  vectors 
are  replaced  by  their  normalized  forms.  Also,  using  (6.42),  (6.43)  and  (6.48),  it  may  be 
established  that  b 2  <  6.  Using  a  similar  approach  as  above,  we  obtain  the  nomalized  split 
lattice  algorithm  based  on  (72),  Such  an  algorithm  consists  in  the  following  computations, 


7*-i 

Ok 


0k  = 

k 

E  = 

q 


N-2 

(  a*^N+k-2,i+k)fik-l 

«=0 

flfc-i  ~  7fc-i 

0  k-  2  —  7*- 2 

2#*-i  -  a*0jfe_2 

*-l  fe-2 

2£-a*£ 

q  q 


(6.75) 
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0k 


h 


fN+k- l{z) 
Pk 


q 

0k 

0k- 1 

1[(1  +  z)fN+k-2(z )  -  Q!fci>fcI1/N+fc-3(^)] 

1  u  o 

1  -  a —  k~2 r-',n 

f'k-l 


The  split  lattice  algorithm  based  on  (72)  requires  2 Nn  +  0.5 n2  +  0(N  +  n)  MULT  and 
3iVn  +  n2  +  0(iV  +  n)  ADD  while  2Nn  +  n2  +  0(N  +  n)  MULT  and  3-/Vn  +  1.5n2  +  0(Af  +  n) 
ADD  are  required  for  the  split  lattice  algorithm  based  on  (6.73).  For  the  normalized  split 
lattice  aigoirthm  described  in  (6.74),  3Arn  +  1.5n2  +  0(A'+n)  MULT,  3iVn  +  1.5n2  +  (N+n) 
ADD  and  n  square-roots  are  required.  The  normalized  split  lattice  algorithm  described  in 
(6.75)  requires  3Nn  +  n2  +  0(N  +  n)  MULT,  3Nn  +  n2  +  0(N  +  n)  ADD  and  n  square-roots. 

The  recurrence  relations  for  the  split  lattice  algorithm  described  here  may  be  employed 
to  obtain  a  new  lattice  realization  of  the  FIR  filter  having  the  transfer  function  Such 

a  realization  is  shown  in  Fig.2  for  order  n=5.  It  is  important  to  note  that  the  previously 
described  split  algorithms  of  [8],  [9]  cannot  be  used  to  realize  the  lattice  implementation 
of  the  FIR  filter  having  the  transfer  function  <j>{z~x),  since  the  terminal  stage  of  such  a 
realization  will  require  division  by  (1  —  z~l),  which  corresponds  to  an  unstable  HR  filter  of 
order  1  having  a  pole  at  z  =  1..  ‘ 
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6.5  Related  Algorithms  and  Extensions 


A  natural  qestion  arises  at  this  stage —  can  similar  savings  in  number  of  arithmetic  op¬ 
erations  be  obtained  for  more  general  forms  of  Toeplitz  matrices,  for  example,  hermitian 
Toeplitz  matrices,  non-symmetric  Toeplitz  matrices,  and  two-dimensional  Toeplitz  matri¬ 
ces?  Also,  another  problem  closely  related  to  the  problem  of  computing  the  Levinson 
polynomial,  is  to  solve  for  the  general  linear  system  Tx  =  b.  In  the  following,  we  briefly 
summarize  some  of  our  findings  on  these  questions.  A  more  detailed  analysis  will  form  a 
topic  of  a  future  research  paper. 

i 

For  a  hermitian  Toeplitz  matrix  T,  the  recurrence  can  be  developed  on  vectors  satisfying 
the  system, 

Tk<ik  =  [Ok  ■  •  •  0k¥  fc  =  l,2,”*,n 

where  Ok  is  a  real  scalar.  Once  again,  q*  satisfies  the  symmetry  property  q*  =  q*,  where 
q*:  is  defined  as  q*  =  «/q£.  In  this  case,  qk,o  ^  1,  k  =  1,2, •  •  •  ,n.  The  recusive  algorithm 
based  on  the  computation  of  q*,  k  =  1, ••*,«,  provides  computational  savings  similar  to 
the  real  symmetric  case  [36]. 

For  a  non-symmetric  Toeplitz  matrix  T,  a  recurrence  may  be  developed  on  vectors  sat¬ 
isfying  the  system  T*q*  =  [Ok  •  •  •  #*]*,  where  Ok  is  a  real  scalar.  However,  in  this  case,  the 
vector  q*  possesses  no  symmetry  properties  and,  as  a  result  ,  no  computational  savings 
are  possible  over  the  Levinson- Durbin  algorithm.  The  derivation  of  such  an  algorithm  is 
straightforward  and  is  omitted  here.  For  the  two-dimensional  case,  there  jure  no  compu¬ 
tational  savings  as  well,  unless  additional  constraints  axe  imposed  on  the  mathematical 
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structure  of  the  matrix. 


The  split  Levinson  algorithm  based  on  recursive  computation  of  symmetric  vectors  q* 
as  described  here  may  also  be  employed  to  solve  for  the  linear  system  Tx  =  b  in  an  efficient 
manner.  An  outline  of  the  algorithm  is  as  follows. 

The  data  vector  b  can  written  as  a  sum  of  a  symmetric  vector  s  and  an  anti-symmetric 
vector  a,  where  s  =  |[b  +  b],  and  a  =  b  —  s  =  |[b  —  b].  If  y  is  solution  to  the  system 
T y  =  s  and  z  is  solution  to  the  system  T z  =  a,  then  solution  to  the  system  Ty  =  b  can 
simply  be  written  as, 

x  =  y  +  z  (6.76) 


Thus,  the  task  of  solving  for  Tx  =  b  is  reduced  to  solving  two  tasks,  that  is,  Ty  =  s 
and  Tz  =  a  at  a  cost  of  2n  ADD  and  0.5n  MULT.  Since  s  and  a  axe  symmetric  and 
antisymmetric  vectors  respectively,  the  corresponding  solution  vectors  y  and  z  are  also 
symmetric  and  antisymmetric  respectively. 

Combining  (6.24)  and  a  scaler  version  of  (6.27)  (  scalar  =  —  gjk-),  and  letting, 


we  get 


Pk  =  qk~  ~~  < 
-2 


Ok 

Ok-2  ^  2 


TkPk  =  fakO  •  •  •  Or*]* 


(6.77) 


(6.78) 


where 


Tk  = 


<Tk 

(1  +  Pk) 


(6.79) 


1 
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Similarly,  recall  from  (6.56)  -  (6.60),  that 


■••(>- rtf 


(6.80) 


where 


P 1  = 


q*-i 

o 


o 

qfc-i 


(6.81) 


and 


<7fc 


(6.82) 


T‘  (i  -  «) 

Given  the  vector  q*,  computation  of  the  vectors  p*  and  pj  using  (6.77)  and  (6.81)  require 
a  total  of  k  MULT  and  k  ADD.  Note  that  the  vectors  p*  and  p£  are  symmetric  and  anti¬ 
symmetric  respectively. 

The  vectors  p*  and  p*  can  be  employed  directly  for  solving  the  systems  Ty  =  s  and 
Tz  =  a  respectively.  The  procedure  for  solving  Ty  =  s  (  or  Tz  =  a  )  consists  in  computing 
the  solution  corresponding  to  all  the  center-subvectors  of  s  (  or  a  ).  Since  all  the  center- 
subvectors  of  s  (  or  a  )  are  symmetric  (  or  anti-symmetric),  the  corresponding  solutions  are 
also  symmetric  (  or  anti-symmetric).  The  algorithms  for  solving  Ty  =  s  using  p*  and  for 
solving  Tz  =  a  using  p£  are  described  in  the  literature  [37],  [38],  and  the  reader  is  referred 
to  them  for  further  details.  Once  z  and  y  are  computed,  the  solution  vector  x  is  obtained 
using  (6.76). 

The  computational  complexity  of  the  algorithm  for  solving  Tx  =  b  using  the  approach 
outlined  above  is  1.125n3-i-0(n)  MULT  and  2n2+0(n)  ADD.  This  represents  a  small  im¬ 
provement  in  computational  complexity  over  the  algorithm  described  in  [39]  which  requires 
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1.25n2  +  0(n )  MULT  and  2n2  -f  0(n )  ADD.  However,  there  is  another  advantage  in  using 
the  approach  described  here.  We  solve  for  the  systems  Ty  =  s  and  T z  =  a  by  recursively 
solving  Tmym  =  sm  and  Tm zm  =  am,  where  sm  and  am  are  center-subvectors  of  s  and  a 
respectively.  Therefore,  ym  and  zm  are  solutions  to  SYM  and  ASYM  part  of  the  center- 
subvectors,  bm,  of  the  data  vector  b  respectively.  Solution  to  the  subsystem  Tmxm  =  bm 
can  simply  be  obtained  as  xm  =  ym  +  zm.  In  contract,  the  algorithm  in  [39]  recursively 
computes  solution  to  the  system  TfeX*,  =  [b^bi  •  •  •  &*]*,  where  ^  b0.  As  a  result,  the 
intermediate  solutions  Xx,  •  •  • ,  x„  are  not  solutions  to  any  subvector  of  data  vector  b.  There¬ 
fore,  the  algorithm  in  [39]  may  not  be  applicable  in  system  applications  that  require  solution 
to  subsystems  as  well. 
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Chapter  7 


Conclusions 


Our  goal  has  been  to  bring  together  the  richness  of  results  available  in  the  theories  of 
orthogonal  polynomials  and  matrices  as  they  relate  to  Levinson,  Schur,  and  lattice  algo¬ 
rithms  and  the  underlying  mathematical  structure  of  the  Toeplitz  matrices.  Emphasis  is 
placed  on  exploiting  the  properties  of  the  Toeplitz  matrices  to  derive  computationally  effi¬ 
cient  algorithms  for  solving  the  associated  linear  system  and  computing  other  parameters 
of  interest. 

In  our  work,  we  have  studied  various  numerical  algorithms  for  Toeplitz  matrices  from  the 
standpoint  of  (i)  computational  complexity,  and  (ii)  numerical  stability  and  finite  precision 
arithmetic  implementation.  In  the  domain  of  computational  complexity,  we  have  derived 
fast  algorithms  for  testing  the  strict  and  wide  sense  stability  of  discrete  time  systems, 
computing  the  reflection  coefficients  and  solving  one  and  two  dimensional  Toeplitz  linear 
system.  A  superfast  algorithm  is  also  derived  for  solving  a  block  Toeplitz  linear  system. 
A  new  lattice  realization  of  a  digital  filter  is  also  obtained.  In  the  domain  of  numerical 
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stability,  it  is  shown  that  the  recently  reported  split  Levinson  algorithm  is  weakly  stable. 
The  finite  precision  implementation  of  the  algorithms  results  in  errors  being  introduced.  We 
have  analyzed  the  Levinson,  Schur  ,  and  lattice  algorithms  for  the  effects  of  finite  precision 
arithmetic  on  their  implementation. 

A  number  of  linear  algebra  problems  related  to  Toeplitz  matrices  remain  open  at  this 
time,  one  of  the  most  significant  being  the  design  of  a  fast  Fourier  transform  based  superfast 
algorithm  for  solving  a  near-to- Toeplitz  linear  system.  Also,  the  application  of  the  new 
algorithms  to  computing  the  eigendecomposition  of  the  Toeplitz  matrices  must  be  examined. 
The  parallel  processing  aspects  of  the  various  algorithms  and  their  suitability  for  VLSI 
implememtation  also  needs  to  studied  for  solving  large  scale  problems. 
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